Pseudomonas Genome Alignment Project

Compare the genomes of Pseudomonas aeruginosa ATCC 27853 and ST235 with quick and dirty but reproducible steps and a dotplot visualization gene.

---

Downloading, aligning, and visualizing the genomes of two Pseudomonas aeruginosa strains using MUMmer, with special focus on identifying and highlighting the VIM gene (blaVIM-2) that is unique to ST235.

---

🛠️ Prerequisites

📦 Install Required Tools

bash sudo apt-get update sudo apt-get install mummer gnuplot ncbi-blast+

1. Download Genome FASTA Files

Download ATCC 27853 genome

ATCC.fasta "https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_001687285.1/"

Download ST235 genome

ST.fasta "https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_016923535.1/"

```bash

Download VIM gene (blaVIM-2) for analysis

curl -L 'https://www.ncbi.nlm.nih.gov/nuccore/GU137304.1?report=fasta&format=text' -o blaVIM2.fasta ```

2. Align Genomes with MUMmer

bash nucmer --prefix=ATCC_vs_ST ATCC.fasta ST.fasta delta-filter -1 ATCC_vs_ST.delta > ATCC_vs_ST.filtered.delta show-coords -rcl ATCC_vs_ST.filtered.delta > ATCC_vs_ST.coords

3. Search for VIM Gene

```bash

Create BLAST databases

makeblastdb -in ATCC.fasta -dbtype nucl makeblastdb -in ST.fasta -dbtype nucl

Search for VIM gene in both genomes

blastn -query blaVIM2.fasta -db ATCC.fasta -out VIMinATCC.txt -outfmt 6 blastn -query blaVIM2.fasta -db ST.fasta -out VIMinST.txt -outfmt 6

Check results

echo "VIM gene hits in ATCC:" wc -l VIMinATCC.txt echo "VIM gene hits in ST235:" wc -l VIMinST.txt ```

4. Summarize Alignment

bash awk 'NR>5 {aligned+=$7; refcov+=$11; querycov+=$12; blocks++} END {print "Total aligned bases:", aligned; print "Alignment blocks:", blocks; print "Avg % identity:", "see below"}' ATCC_vs_ST.coords awk 'NR>5 {id+=$8*$7; len+=$7} END {if(len>0) print "Weighted avg % identity:", id/len; else print "No alignments"}' ATCC_vs_ST.coords

5. Generate Basic Dotplot

```bash mummerplot --png --large --layout --filter --prefix=ATCCvsSTplot ATCCvs_ST.filtered.delta

Fix gnuplot error if it occurs

sed -i '/set mouse clipboardformat/d' ATCCvsSTplot.gp gnuplot ATCCvsSTplot.gp ```

6. Create Enhanced Dotplot with VIM Gene Highlight

```bash

Add VIM gene annotation and styling to the gnuplot script

cat >> ATCCvsST_plot.gp << 'EOF'

Clear legend and annotations

set key at graph 0.02, graph 0.98 left font "Arial,12"

Add red arrow pointing to VIM gene location

set arrow 1 from graph 0.9, first 7074888 to graph 0.95, first 7074888 lc rgb "red" lw 3 head filled set label 1 "VIM gene" at graph 0.85, first 7074888 center tc rgb "red" font "Arial,12,bold"

Add grey dots for unique regions

set object 2 circle at graph 0.99, first 7072697 size graph 0.005 fc rgb "#808080" fillstyle solid noborder set object 3 circle at first 50000, graph 0.005 size graph 0.005 fc rgb "#808080" fillstyle solid noborder set object 4 circle at first 650000, graph 0.005 size graph 0.005 fc rgb "#808080" fillstyle solid noborder set object 5 circle at first 1250000, graph 0.005 size graph 0.005 fc rgb "#808080" fillstyle solid noborder

Add title and better axis labels

set title "Genome Alignment: ATCC 27853 vs ST235\nCyan=Forward alignments, Purple=Reverse alignments, Red=VIM gene, Grey=Unique regions" font "Arial,14,bold" set xlabel "ATCC 27853 genome position" font "Arial,12" set ylabel "ST235 genome position" font "Arial,12" EOF

Regenerate the enhanced plot

gnuplot ATCCvsST_plot.gp ```

Results Summary

| Metric | Value | |-----------------------------|--------------| | Total aligned bases | 6,123,079 | | Alignment blocks | 298 | | Weighted avg % identity | >98% | | VIM gene in ATCC | Not found | | VIM gene in ST235 | Found (contig JAFFXY010000040.1) |

Key Findings:

🖼️ Enhanced Alignment Dotplot

Plot Legend:

Plot legend (colors, meaning and counts)

• FWD (cyan): forward alignments; here 169 blocks covering 3,470,819 bp of the reference.
• REV (purple): reverse/inverted alignments; here 129 blocks covering 2,652,260 bp of the reference.
• OTH (black): other/auxiliary alignments or overlays (used as a placeholder in the script); 0 blocks in this run.

Notes: - The plot line-styles have been explicitly set in pa_comparison.gp so plotted colors match the legend (ls 1 = cyan, ls 2 = purple, ls 3 = black). - "OTH" is a generic label for any additional alignment type or overlay plotted with the third line style; it is present as a placeholder so the legend can be extended without reformatting the plot. - The visual prominence of purple dots on the dotplot can be affected by placement (reverse matches are often scattered and visually distinct), overplotting, and color contrast — the numeric counts above show forward alignments still cover more bases in this comparison.

📚 References & Further Info

🧑‍🔬 Legacy/Alternative Data & Methods

Status

• ATCC 27853 complete genome: CP015117.1 → ATCC_27853_CP015117.1.fasta (NCBI)
• ATCC 27853 alternative genome: CP011857.1 → ATCC_27853_alternative.fasta (NCBI) - from PMC5467263
• ST235 representative complete genome (NCGM2.S1): AP012280.1 → ST235_NCGM2S1_AP012280.1.fasta (NCBI)
• Note: Urbanowicz et al., 2021 link LFMO00000000 is a WGS master and contains no sequence data (NCBI) ```bash # Download ATCC 27853 genome (CP015117.1) wget -O ATCC27853CP015117.1.fasta "https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=nuccore&report=fasta&id=CP015117.1"

# Download ATCC 27853 alternative genome (CP011857.1) from PMC5467263 wget -O ATCC27853alternative.fasta "https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=nuccore&report=fasta&id=CP011857.1"

# Download ST235 genome (NCGM2.S1) wget -O ST235NCGM2S1AP012280.1.fasta "https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=nuccore&report=fasta&id=AP012280.1"

# Download and extract minimap2 wget -O minimap2-2.28x64-linux.tar.bz2 https://github.com/lh3/minimap2/releases/download/v2.28/minimap2-2.28x64-linux.tar.bz2 tar -xjf minimap2-2.28_x64-linux.tar.bz2

# Run genome alignment (CP015117.1 vs ST235) ./minimap2-2.28x64-linux/minimap2 -x asm5 -t 4 ATCC27853CP015117.1.fasta ST235NCGM2S1AP012280.1.fasta > ATCC27853vsST235NCGM2S1.paf

# Run genome alignment (CP011857.1 vs ST235) ./minimap2-2.28x64-linux/minimap2 -x asm5 -t 4 ATCC27853alternative.fasta ST235NCGM2S1AP012280.1.fasta > ATCC27853CP011857vsST235_NCGM2S1.paf

# Quick summary statistics awk 'BEGIN{m=0;b=0;ql=0;tl=0} {m+=$10; b+=$11; if($2>ql) ql=$2; if($7>tl) tl=$7} END{printf "matches=%d alnbases=%d qlen=%d tlen=%d pid=%.6f\n",m,b,ql,tl,m/b}' ATCC27853vsST235NCGM2S1.paf awk 'BEGIN{m=0;b=0;ql=0;tl=0} {m+=$10; b+=$11; if($2>ql) ql=$2; if($7>tl) tl=$7} END{printf "matches=%d alnbases=%d qlen=%d tlen=%d pid=%.6f\n",m,b,ql,tl,m/b}' ATCC27853CP011857vsST235NCGM2S1.paf `` - Tool:minimap2(v2.28) with-x asm5 - Output files: -ATCC27853vsST235NCGM2S1.paf(CP015117.1 vs ST235) -ATCC27853CP011857vsST235NCGM2S1.paf(CP011857.1 vs ST235) - Quick summary (PAF aggregate): - CP015117.1 vs ST235: matches = 5,143,335; aligned bases = 7,093,386; rough pid = 0.725 - CP011857.1 vs ST235: matches = 5,141,970; aligned bases = 7,146,480; rough pid = 0.720 - If available, provide the exact ST235 accession from Urbanowicz et al., 2021 to replace the representative genome and re-run. - Compute robust ANI and SNP/indel stats (e.g.,dnadiff/fastANI`) and generate a brief report.

### Plot sanitization (headless rendering)

Some gnuplot scripts produced by the plot generator include interactive commands (mouse/clipboard handlers, print banners and pause -1) which force gnuplot into interactive mode and prevent headless PNG generation. To run the plots on a headless machine or inside scripts, the following conservative changes were applied to problematic .gp files in this repo:

• Removed interactive lines:
  • any lines beginning with set mouse ...
  • print lines that show an interactive banner
  • pause -1
• Ensured non-interactive output by adding (when needed):
  • set terminal pngcairo size 1600,1200 enhanced font "Courier,8"
  • set output '/full/path/to/<name>.png' (absolute path) before plotting commands

Files already sanitized in this repo: - pa_comparison.gp → pa_comparison.png - pa_comparison.relaxed.gp → pa_comparison.relaxed.png

How to re-run the sanitized plots - From the project root run:

bash cd /home/ash022/pseudomonas_align_strains gnuplot pa_comparison.gp gnuplot pa_comparison.relaxed.gp gnuplot pa_comparison.rects.gp

Quick check (verify output):

bash ls -l pa_comparison*.png # Expect non-zero-sized PNG files: pa_comparison.png, pa_comparison.relaxed.png, pa_comparison.rects.png

Run the rect generator (Makefile & wrapper)

There are convenient, safe ways to generate the rectangle overlays and vector exports without accidentally overwriting a hand-edited pa_comparison.rects.gp.

From the project root you can run the Makefile targets:

```bash

Safe (writes auto-named files, won't overwrite a manual pa_comparison.rects.gp):

make rects

Force overwrite canonical files (only use if you want to replace a hand-edited pa_comparison.rects.gp):

make rects-force ```

Or use the wrapper script which also writes auto-named outputs:

bash scripts/run_make_rects.sh

If you prefer to run the script directly and control filenames/flags, here's the full CLI example that produces PNG, SVG and PDF:

bash python3 scripts/make_rects.py --out-gp pa_comparison.rects.gp \ --out-png pa_comparison.rects.png --export-svg --export-pdf

Note: the default (safe) targets create pa_comparison.rects.auto.gp / .auto.png / .auto.svg / .auto.pdf so you can inspect the generated plots before overwriting any canonical files.

Optional: small helper to sanitize other .gp files - Save this as scripts/sanitize_gp.sh, make it executable, then run it on any .gp files you need to sanitize.

```bash mkdir -p scripts cat > scripts/sanitize_gp.sh <<'EOF' #!/usr/bin/env bash set -euo pipefail for f in "$@"; do echo "Sanitizing: $f" # remove interactive lines that break headless rendering sed -i '/^set mouse/d;/^print /d;/^pause -1/d' "$f"

# ensure a pngcairo terminal is present
if ! grep -q '^set terminal' "$f"; then
  sed -i '1i set terminal pngcairo size 1600,1200 enhanced font "Courier,8"' "$f"
fi

# ensure an absolute output line is present
if ! grep -q '^set output' "$f"; then
  out="$(pwd)/${f%.gp}.png"
  sed -i "1i set output '$out'" "$f"
fi

done EOF chmod +x scripts/sanitize_gp.sh

# Usage example: # scripts/sanitizegp.sh pacomparison.gp pacomparison.relaxed.gp otherplot.gp ```

Export high-resolution figures (SVG / PDF)

The plots generated by mummerplot/gnuplot can be exported as publication-ready SVG or PDF files by adjusting the gnuplot terminal and re-running the script. Example steps I used in this project:

1. Make an SVG copy of a generated .gp and set a high-resolution terminal and output:

bash cp pa_comparison.relaxed.gp pa_comparison.relaxed.svg.gp sed -i "s#^set terminal.*#set terminal svg size 2000,1200#" pa_comparison.relaxed.svg.gp sed -i "s#^set output.*#set output 'pa_comparison.relaxed.svg'#" pa_comparison.relaxed.svg.gp gnuplot pa_comparison.relaxed.svg.gp

1. Make a PDF copy using pdfcairo:

bash cp pa_comparison.relaxed.gp pa_comparison.relaxed.pdf.gp sed -i "s#^set terminal.*#set terminal pdfcairo size 11in,8in enhanced font 'Arial,12'#" pa_comparison.relaxed.pdf.gp sed -i "s#^set output.*#set output 'pa_comparison.relaxed.pdf'#" pa_comparison.relaxed.pdf.gp gnuplot pa_comparison.relaxed.pdf.gp

Files produced in this run (added to the repository):

• pa_comparison.relaxed.svg — high-resolution SVG of the relaxed plot
• pa_comparison.relaxed.pdf — PDF export of the relaxed plot
• pa_comparison.rects.svg / pa_comparison.rects.pdf — precise-rectangles plot (SVG + PDF)
• pa_comparison.annotated3.svg / pa_comparison.annotated3.pdf — annotated version (SVG + PDF)

Notes: - These copies preserve the original .gp script; we create .svg.gp / .pdf.gp variants so the original scripts remain unchanged. - SVG is preferred for vector editing and publication; PDF via pdfcairo is suitable for inclusion in manuscripts.

Reproducible pipeline — commands and helper scripts

Below are concrete commands and the small helper scripts used during this analysis so you can reproduce the pipeline from Bakta TSVs to grouped summaries and a proportional Venn diagram.

Prereqs (system packages) ```bash

required system tools used in this repo

sudo apt-get update sudo apt-get install -y mummer gnuplot ncbi-blast+ python3-pip ```

Prereqs (python packages) ```bash

install plotting and venn library into the same python used here

python3 -m pip install --user matplotlib matplotlib-venn ```

Assumptions - Input annotation TSVs are ATCC.tsv and ST.tsv (Bakta output TSVs). If you don't have them, run Bakta or Prokka first to create those files.

1) Group each Bakta TSV by the raw Product string and add a Matches count

Save this as scripts/group_by_product_with_counts.py and run it.

```python

!/usr/bin/env python3

from pathlib import Path from collections import defaultdict

def groupbyproduct(inpath, outpath): p = Path(inpath) text = p.readtext(encoding='utf-8', errors='surrogateescape').splitlines() headerline = None for i, l in enumerate(text[:10]): if l.startswith('#') and 'Sequence' in l: headerline = l.lstrip('#') start = i break if headerline is None: headerline = text[0] start = 0 fieldnames = [f.strip() for f in headerline.split('\t')] rows = [] for l in text[start+1:]: if not l.strip(): continue parts = l.split('\t') if len(parts) < len(fieldnames): parts += [''] * (len(fieldnames) - len(parts)) row = dict(zip(fieldnames, parts)) rows.append(row) groups = defaultdict(list) for r in rows: prod = r.get('Product','') groups[prod].append(r) outheader = ['Product','Matches'] + [c for c in fieldnames if c!='Product'] outlines = [] for prod, items in sorted(groups.items(), key=lambda x: (-len(x[1]), x[0])): matches = len(items) agg = {} for col in fieldnames: vals = [] for it in items: v = it.get(col,'') if v and v not in vals: vals.append(v) agg[col] = ';'.join(vals) rowout = [prod, str(matches)] + [agg[c] for c in fieldnames if c!='Product'] outlines.append('\t'.join(rowout)) outp = Path(outpath) outp.parent.mkdir(parents=True, existok=True) with outp.open('w', encoding='utf-8') as fh: fh.write('\t'.join(outheader) + '\n') for l in outlines: fh.write(l + '\n')

if name == 'main': groupbyproduct('ATCC.tsv','ATCCgroupedbyproductwithcounts.tsv') groupbyproduct('ST.tsv','STgroupedbyproductwithcounts.tsv') ```

Run: bash python3 scripts/group_by_product_with_counts.py

2) Merge the two grouped files keeping per-file counts and a Total_Matches column

Save this as scripts/merge_grouped_with_counts.py.

```python

!/usr/bin/env python3

from pathlib import Path

def readgrouped(path): text = Path(path).readtext(encoding='utf-8', errors='surrogateescape').splitlines() header = text[0].split('\t') rows = {} for ln in text[1:]: if not ln.strip(): continue parts = ln.split('\t') if len(parts) < len(header): parts += [''] * (len(header)-len(parts)) d = dict(zip(header, parts)) prod = d.get('Product','') rows[prod] = {'Matches': int(d.get('Matches','0') or 0), 'data': d} return header, rows

ha, ra = readgrouped('ATCCgroupedbyproductwithcounts.tsv') hs, rs = readgrouped('STgroupedbyproductwithcounts.tsv') colsa = [c for c in ha if c not in ('Product','Matches')] colss = [c for c in hs if c not in ('Product','Matches')] allproducts = sorted(set(list(ra.keys()) + list(rs.keys()))) outcols = ['Product','ATCCMatches','STMatches','TotalMatches'] + [f'ATCC{c}' for c in colsa] + [f'ST{c}' for c in colss] with open('combinedgroupedwithcounts.tsv','w',encoding='utf-8') as fh: fh.write('\t'.join(outcols) + '\n') for prod in allproducts: a = ra.get(prod) s = rs.get(prod) amatches = a['Matches'] if a else 0 smatches = s['Matches'] if s else 0 total = amatches + smatches row = [prod, str(amatches), str(smatches), str(total)] for c in colsa: row.append(a['data'].get(c,'') if a else '') for c in colss: row.append(s['data'].get(c,'') if s else '') fh.write('\t'.join(row) + '\n') ```

Run: bash python3 scripts/merge_grouped_with_counts.py

3) Mark source (ATCConly / STonly / both)

Save as scripts/mark_sources.py:

```python

!/usr/bin/env python3

from pathlib import Path lines = Path('combinedgroupedwithcounts.tsv').readtext(encoding='utf-8', errors='surrogateescape').splitlines() header = lines[0].split('\t') with open('combinedgroupedwithcountsmarked.tsv','w',encoding='utf-8') as fh: fh.write('\t'.join(header + ['Source']) + '\n') for ln in lines[1:]: if not ln.strip(): continue parts = ln.split('\t') row = dict(zip(header, parts)) a = int(row.get('ATCCMatches','0') or 0) s = int(row.get('STMatches','0') or 0) if a>0 and s>0: src='both' elif a>0: src='ATCConly' elif s>0: src='STonly' else: src='none' fh.write('\t'.join(parts + [src]) + '\n') ```

Run: bash python3 scripts/mark_sources.py

4) Proportional two-set Venn (matplotlib-venn)

Save as scripts/venn_proportional.py:

```python

!/usr/bin/env python3

from pathlib import Path from matplotlib_venn import venn2 import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt

lines = Path('combinedgroupedwithcountsmarked.tsv').readtext(encoding='utf-8', errors='surrogateescape').splitlines() header = lines[0].split('\t') srcidx = header.index('Source') counts = {'ATCConly':0,'STonly':0,'both':0} for ln in lines[1:]: if not ln.strip(): continue parts = ln.split('\t') src = parts[srcidx] if src in counts: counts[src]+=1 set1 = counts['ATCConly'] + counts['both'] set2 = counts['STonly'] + counts['both'] inter = counts['both'] plt.figure(figsize=(6,6)) venn = venn2(subsets=(set1-inter, set2-inter, inter), setlabels=('ATCC','ST')) if venn.getlabelbyid('10'): venn.getlabelbyid('10').settext(str(counts['ATCConly'])) if venn.getlabelbyid('01'): venn.getlabelbyid('01').settext(str(counts['STonly'])) if venn.getlabelbyid('11'): venn.getlabelbyid('11').settext(str(counts['both'])) plt.title('Product overlap: ATCC vs ST') plt.savefig('vennproportionalATCCST.png', dpi=200, bbox_inches='tight') ```

Run: bash python3 scripts/venn_proportional.py

Files created by the reproducible pipeline (names used above) - ATCC_grouped_by_product_with_counts.tsv - ST_grouped_by_product_with_counts.tsv - combined_grouped_with_counts.tsv - combined_grouped_with_counts_marked.tsv - venn_proportional_ATCC_ST.png (proportional Venn diagram)

If you want, I can add and commit these helper scripts into the repository under scripts/ now. Reply "yes, commit scripts" to have me create the scripts/ folder and add these files to git.

Bakta installation & annotation (how I ran it)

The genome annotation step in this project used Bakta. Below are the concrete commands I used so you can reproduce the annotation locally.

Prereqs (recommended) ```bash

optional: create an isolated conda env (recommended)

conda create -n bakta python=3.10 -y conda activate bakta

or use your system Python/pip

python3 -m pip install --user pipx || true ```

Install Bakta ```bash

using pip inside the active environment

python3 -m pip install --upgrade pip python3 -m pip install bakta

confirm version

bakta --version

expected output (example): bakta 1.11.3

```

Obtain a Bakta DB bundle

• I used a lightweight DB bundle (named db-light.tar.xz in my run). If you have a supplied archive use it; otherwise download the appropriate Bakta DB from the Bakta release or provider.

```bash

create a folder for the DB and extract the bundle

mkdir -p bakta_db

path/to/db-light.tar.xz may be a local path or downloaded file

tar -xJf /path/to/db-light.tar.xz -C bakta_db ```

Run Bakta annotation (examples used in this repository) ```bash

annotate ATCC

bakta annotate \ --db-dir baktadb \ --output baktaATCC \ --prefix ATCC \ ATCC.fasta

annotate ST

bakta annotate \ --db-dir baktadb \ --output baktaST \ --prefix ST \ ST.fasta ```

Outputs - Annotation directories: bakta_ATCC/, bakta_ST/ - Main TSVs used downstream: ATCC.tsv and ST.tsv (extracted from the Bakta outputs)

Troubleshooting notes - If Bakta (or one of its bundled AMR binaries) fails with C++ runtime errors (missing GLIBCXX symbols), use a Conda environment with a current toolchain or update the environment's libstdc++/gcc packages. For example: bash conda install -n bakta libgcc-ng libstdcxx-ng -y - If you have a custom DB bundle, pass its directory with --db-dir as shown above.

If you want, I can add the exact bakta command-lines that were run (including any extra flags) into the repository's run-log; tell me if you'd like them appended.

Exact commands I ran for Bakta in this project

The following are the exact shell commands used during this session to set up Bakta, install the DB bundle provided to me, fix the AMR runtime libraries, run annotation, and copy the produced TSVs into the repository root.

```bash

create & enter conda env (optional but recommended)

conda create -n bakta python=3.10 -y conda activate bakta

install bakta

python3 -m pip install --upgrade pip python3 -m pip install bakta

extract the user-supplied Bakta DB bundle used in this run

mkdir -p bakta_db

(path used in this session)

replace the path below with your DB archive path if different

tar -xJf /mnt/z/Download/db-light.tar.xz -C bakta_db

fix runtime libs for bundled AMR binary (if you see GLIBCXX errors)

conda install -n bakta libgcc-ng libstdcxx-ng -y

Run Bakta annotation (commands used here)

bakta annotate --db-dir baktadb --output baktaATCC --prefix ATCC ATCC.fasta bakta annotate --db-dir baktadb --output baktaST --prefix ST ST.fasta

Copy or move main TSV outputs into repo root for downstream scripts

(Bakta writes TSVs into its output directory; copy them here)

cp baktaATCC/ATCC.tsv ./ATCC.tsv cp baktaST/ST.tsv ./ST.tsv ```

These exact commands should reproduce the annotation step performed for ATCC.fasta and ST.fasta in this project. If your DB archive or paths differ, update the tar -xJf command and --db-dir accordingly.

Imported Bakta TSVs, grouping and Venn diagram

I copied Bakta TSV exports (baktaATCC.tsv, baktaST.tsv) into the repository root and produced grouped summaries and a proportional Venn diagram for ATCC vs ST products.

Files created in this step (present in the repo):

• baktaATCC.tsv — Bakta TSV for ATCC
• baktaST.tsv — Combined Bakta TSV(s) for ST
• baktaATCC_grouped_by_product_with_counts.tsv
• baktaST_grouped_by_product_with_counts.tsv
• combined_grouped_with_counts.tsv
• combined_grouped_with_counts_marked.tsv (has Source column: ATCConly/STonly/both)
• venn_proportional_ATCC_ST.png

Reproduce these steps locally (from project root):

```bash

Ensure you have the Bakta TSV exports in the repo root as baktaATCC.tsv and baktaST.tsv.

Group by product (the helper groups by the raw Product string):

python3 - <<'PY' from scripts.groupbyproductwithcounts import groupbyproduct groupbyproduct('baktaATCC.tsv','baktaATCCgroupedbyproductwithcounts.tsv') groupbyproduct('baktaST.tsv','baktaSTgroupedbyproductwithcounts.tsv') PY

Copy grouped files to the canonical names expected by the merge script and merge:

cp baktaATCCgroupedbyproductwithcounts.tsv ATCCgroupedbyproductwithcounts.tsv cp baktaSTgroupedbyproductwithcounts.tsv STgroupedbyproductwithcounts.tsv python3 scripts/mergegroupedwith_counts.py

Mark sources (adds 'Source' column)

python3 scripts/mark_sources.py

Create proportional Venn diagram (requires matplotlib-venn)

python3 scripts/venn_proportional.py ```

Notes:

• Install matplotlib-venn if missing: python3 -m pip install --user matplotlib-venn.
• These helper scripts are lightweight; you may want to run normalization (lowercase, synonym mapping) on Product before grouping for cleaner results.

NOTE: these bakta*.tsv files were produced by the online Proksee/Bakta service (exported from Proksee projects):

• ST combined from Proksee projects:
  • https://proksee.ca/projects/60766a46-90fd-43ae-acce-c5f56372087a
  • https://proksee.ca/projects/cdf82d4c-7de5-4e96-afae-9812033db694
• ATCC from Proksee project:
  • https://proksee.ca/projects/f60d33bc-d8d6-478e-a2a9-93251a45dc59

If you'd like to compare these online-generated Bakta TSVs with local Bakta results, run Bakta locally and recreate the TSVs (see the "Bakta installation & annotation" section above), then re-run the grouping/merge/venn steps and diff the outputs. Example quick comparison commands:

```bash

run local Bakta (example)

bakta annotate --db-dir baktadb --output baktalocalATCC --prefix ATCC ATCC.fasta cp baktalocalATCC/ATCC.tsv ./baktalocal_ATCC.tsv

produce grouped file locally

python3 - <<'PY' from scripts.groupbyproductwithcounts import groupbyproduct groupbyproduct('baktalocalATCC.tsv','baktalocalATCCgroupedbyproductwithcounts.tsv') PY

quick diff against the Proksee-exported grouped file

diff -u baktaATCCgroupedbyproductwithcounts.tsv baktalocalATCCgroupedbyproductwith_counts.tsv | sed -n '1,200p' ```

Comparing online vs local Bakta runs is useful to detect differences in DB versions, annotation policies, or minor version-dependent behavior in AMR/feature detection.

Genome alignment with NUCmer and visualization

The repository contains an example whole-genome alignment workflow using MUMmer's nucmer plus downstream utilities (delta-filter, show-coords, show-snps) and dnadiff for a summary. The generated pa_comparison.* files in the project are the artifacts produced by that run. The minimal command sequence to reproduce the alignment and produce dotplots is shown below (run from the project root):

```bash

Run nucmer (ATCC as reference, ST as query)

nucmer --prefix=pa_comparison ATCC.fixed.fasta ST.fixed.fasta

Optional: filter delta to 1-to-1 best matches

delta-filter -1 pacomparison.delta > pacomparison.filtered.delta || true

Generate alignment coordinates

show-coords -rcl pacomparison.filtered.delta > pacomparison.coords

SNP/indel summary

show-snps -Clr pacomparison.filtered.delta > pacomparison.snps

Whole-genome summary using dnadiff (optional)

dnadiff -p pa_comparison ATCC.fixed.fasta ST.fixed.fasta ```

Visualization (dotplot using mummerplot and gnuplot):

```bash

Prefer the filtered delta file when available; fall back to the unfiltered delta.

if [ -s pacomparison.filtered.delta ]; then mummerplot --fat --layout --filter -p pacomparison pacomparison.filtered.delta || mummerplot --fat --layout -p pacomparison pacomparison.filtered.delta elif [ -s pacomparison.delta ]; then mummerplot --fat --layout -p pacomparison pacomparison.delta || true else echo "No delta file found: pacomparison.filtered.delta or pacomparison.delta" >&2 fi

Run the generated gnuplot script in headless mode (the repo contains sanitized .gp scripts)

gnuplot pa_comparison.gp ```

Notes: - The repo already contains sanitized and exported dotplot files (pa_comparison.png, pa_comparison.svg, pa_comparison.pdf, pa_comparison.rects.* etc.). - If mummerplot or the MUMmer tools are not on your PATH, install MUMmer (the version used here included nucmer, delta-filter, show-coords, show-snps, mummerplot, and dnadiff).

I will commit the pa_comparison.* artifacts that exist in the repository so the outputs are tracked. If you prefer a smaller set (for example only pa_comparison.coords, pa_comparison.png, pa_comparison.pdf, and pa_comparison.snps) tell me and I'll commit only those.

Results — comparison artifacts

The repository now contains a small set of comparison artifacts created from two snapshots of combined_grouped_with_counts_marked.tsv (the original and the latest committed versions). These files are stored under results/ for provenance and easy inspection:

• results/orig_snapshot.tsv — the original combined_grouped_with_counts_marked.tsv (first committed version).
• results/latest_snapshot.tsv — the latest combined_grouped_with_counts_marked.tsv (HEAD at the time of comparison).
• results/added_products.csv — products present in the latest snapshot but not in the original; includes columns: Product, ATCCMatches, STMatches, Total_Matches, Source.
• results/removed_products.csv — products present in the original snapshot but not in the latest; same columns.

These CSVs were generated to make it easy to review annotation changes (newly-added products and products that were removed between the two snapshots). If you want to reproduce the comparison locally, one simple way is to extract the two snapshots from Git and re-run the text-based diff pipeline. For example:

```bash

extract two committed versions of the file (replace and with the desired refs)

then run a minimal set-difference on the Product column

git show :combinedgroupedwithcountsmarked.tsv > /tmp/origsnapshot.tsv git show :combinedgroupedwithcountsmarked.tsv > /tmp/latestsnapshot.tsv cut -f1 /tmp/origsnapshot.tsv | tail -n +2 | sort > /tmp/origproducts.txt cut -f1 /tmp/latestsnapshot.tsv | tail -n +2 | sort > /tmp/latestproducts.txt comm -23 /tmp/latestproducts.txt /tmp/origproducts.txt > /tmp/addedproducts.txt comm -13 /tmp/latestproducts.txt /tmp/origproducts.txt > /tmp/removedproducts.txt ```

If you'd like, I can also add the exact commands and commit hashes used to produce the results/ files in this project (and commit them) so the extraction above is completely reproducible; tell me if you want me to add that detail to this section.