1. Installation and Database Setup
To install and prepare NCBI AMRFinderPlus in the bacto environment:
mamba activate bacto
mamba install ncbi-amrfinderplus
mamba update ncbi-amrfinderplus
mamba activate bacto
amrfinder -u
- This will:
- Download and install the latest AMRFinderPlus version and its database.
- Create
/home/jhuang/mambaforge/envs/bacto/share/amrfinderplus/data/.
- Symlink the latest database version for use.
Check available organism options for annotation:
amrfinder --list_organisms
- Supported values include species such as
Escherichia, Klebsiella_pneumoniae, Enterobacter_cloacae, Pseudomonas_aeruginosa and many others.
2. Batch Analysis: Bash Script for Genome Screening
Use the following script to screen multiple genomes using AMRFinderPlus and output only β-lactam/beta-lactamase hits from a metadata table.
Input:
genome_metadata.tsv — tab-separated columns: filename_TAB_organism, with header.
Run:
cd ~/DATA/Data_Patricia_AMRFinderPlus_2025/genomes
./run_amrfinder_beta_lactam.sh genome_metadata.tsv
Script logic:
- Validates metadata input and AMRFinder installation.
- Loops through each genome in the metadata table:
- Maps text organism names to proper AMRFinder
--organism codes when possible (“Escherichia coli” → --organism Escherichia).
- Executes AMRFinderPlus, saving output for each isolate.
- Collects all individual output tables.
- After annotation, Python code merges all results, filters for β-lactam/beta-lactamase genes, and creates summary tables.
Script:
#!/usr/bin/env bash
set -euo pipefail
META_FILE="${1:-}"
if [[ -z "$META_FILE" || ! -f "$META_FILE" ]]; then
echo "Usage: $0 genome_metadata.tsv" >&2
exit 1
fi
OUTDIR="amrfinder_results"
mkdir -p "$OUTDIR"
echo ">>> Checking AMRFinder installation..."
amrfinder -V || { echo "ERROR: amrfinder not working"; exit 1; }
echo
echo ">>> Running AMRFinderPlus on all genomes listed in $META_FILE"
# --- loop over metadata file ---
# expected columns: filename
<TAB>organism
tail -n +2 "$META_FILE" | while IFS=$'\t' read -r fasta organism; do
# skip empty lines
[[ -z "$fasta" ]] && continue
if [[ ! -f "$fasta" ]]; then
echo "WARN: FASTA file '$fasta' not found, skipping."
continue
fi
isolate_id="${fasta%.fasta}"
# map free-text organism to AMRFinder --organism names (optional)
org_opt=""
case "$organism" in
"Escherichia coli") org_opt="--organism Escherichia" ;;
"Klebsiella pneumoniae") org_opt="--organism Klebsiella_pneumoniae" ;;
"Enterobacter cloacae complex") org_opt="--organism Enterobacter_cloacae" ;;
"Citrobacter freundii") org_opt="--organism Citrobacter_freundii" ;;
"Citrobacter braakii") org_opt="--organism Citrobacter_freundii" ;;
"Pseudomonas aeruginosa") org_opt="--organism Pseudomonas_aeruginosa" ;;
# others (Providencia stuartii, Klebsiella aerogenes)
# currently have no organism-specific rules in AMRFinder, so we omit --organism
*) org_opt="" ;;
esac
out_tsv="${OUTDIR}/${isolate_id}.amrfinder.tsv"
echo " - ${fasta} (${organism}) -> ${out_tsv} ${org_opt}"
amrfinder -n "$fasta" -o "$out_tsv" --plus $org_opt
done
echo ">>> AMRFinderPlus runs finished. Filtering β-lactam hits..."
python3 - "$OUTDIR" << 'EOF'
import sys, os, glob
outdir = sys.argv[1]
files = sorted(glob.glob(os.path.join(outdir, "*.amrfinder.tsv")))
if not files:
print("ERROR: No AMRFinder output files found in", outdir)
sys.exit(1)
try:
import pandas as pd
use_pandas = True
except ImportError:
use_pandas = False
def read_one(path):
import pandas as _pd
# AMRFinder TSV is tab-separated with a header line
df = _pd.read_csv(path, sep='\t', dtype=str)
df.columns = [c.strip() for c in df.columns]
# add isolate_id from filename
isolate_id = os.path.basename(path).replace(".amrfinder.tsv", "")
df["isolate_id"] = isolate_id
return df
if not use_pandas:
print("WARNING: pandas not installed; only raw TSV merging will be done.")
# very minimal merging: just concatenate files
with open("beta_lactam_all.tsv", "w") as out:
first = True
for f in files:
with open(f) as fh:
header = fh.readline()
if first:
out.write(header.strip() + "\tisolate_id\n")
first = False
for line in fh:
if not line.strip():
continue
iso = os.path.basename(f).replace(".amrfinder.tsv", "")
out.write(line.rstrip("\n") + "\t" + iso + "\n")
sys.exit(0)
# --- full pandas-based processing ---
dfs = [read_one(f) for f in files]
df = pd.concat(dfs, ignore_index=True)
# normalize column names (lowercase, no spaces) for internal use
norm_cols = {c: c.strip().lower().replace(" ", "_") for c in df.columns}
df.rename(columns=norm_cols, inplace=True)
# try to locate key columns with flexible names
def pick(*candidates):
for c in candidates:
if c in df.columns:
return c
return None
col_gene = pick("gene_symbol", "genesymbol")
col_seq = pick("sequence_name", "sequencename")
col_class = pick("class")
col_subcls = pick("subclass")
col_ident = pick("%identity_to_reference_sequence", "identity")
col_cov = pick("%coverage_of_reference_sequence", "coverage_of_reference_sequence")
col_iso = "isolate_id"
missing = [c for c in [col_gene, col_seq, col_class, col_subcls, col_iso] if c is None]
if missing:
print("ERROR: Some required columns are missing in AMRFinder output:", missing)
sys.exit(1)
# β-lactam filter: class==AMR and subclass contains "beta-lactam"
mask = (df[col_class].str.contains("AMR", case=False, na=False) &
df[col_subcls].str.contains("beta-lactam", case=False, na=False))
df_beta = df.loc[mask].copy()
if df_beta.empty:
print("WARNING: No β-lactam hits found.")
else:
print(f"Found {len(df_beta)} β-lactam / β-lactamase hits.")
# write full β-lactam table
beta_all_tsv = "beta_lactam_all.tsv"
df_beta.to_csv(beta_all_tsv, sep='\t', index=False)
print(f">>> β-lactam / β-lactamase hits written to: {beta_all_tsv}")
# -------- summary by gene (with list of isolates) --------
group_cols = [col_gene, col_seq, col_subcls]
def join_isolates(vals):
# unique, sorted isolates as comma-separated string
uniq = sorted(set(vals))
return ",".join(uniq)
summary_gene = (
df_beta
.groupby(group_cols, dropna=False)
.agg(
n_isolates=(col_iso, "nunique"),
isolates=(col_iso, join_isolates),
n_hits=("isolate_id", "size")
)
.reset_index()
)
# nicer column names for output
summary_gene.rename(columns={
col_gene: "Gene_symbol",
col_seq: "Sequence_name",
col_subcls: "Subclass"
}, inplace=True)
sum_gene_tsv = "beta_lactam_summary_by_gene.tsv"
summary_gene.to_csv(sum_gene_tsv, sep='\t', index=False)
print(f">>> Gene-level summary written to: {sum_gene_tsv}")
print(" (now includes 'isolates' = comma-separated isolate_ids)")
# -------- summary by isolate & gene (with annotation) --------
agg_dict = {
col_gene: ("Gene_symbol", "first"),
col_seq: ("Sequence_name", "first"),
col_subcls: ("Subclass", "first"),
}
if col_ident:
agg_dict[col_ident] = ("%identity_min", "min")
agg_dict[col_ident + "_max"] = ("%identity_max", "max")
if col_cov:
agg_dict[col_cov] = ("%coverage_min", "min")
agg_dict[col_cov + "_max"] = ("%coverage_max", "max")
# build aggregation manually (because we want nice column names)
gb = df_beta.groupby([col_iso, col_gene], dropna=False)
rows = []
for (iso, gene), sub in gb:
row = {
"isolate_id": iso,
"Gene_symbol": sub[col_gene].iloc[0],
"Sequence_name": sub[col_seq].iloc[0],
"Subclass": sub[col_subcls].iloc[0],
"n_hits": len(sub)
}
if col_ident:
vals = pd.to_numeric(sub[col_ident], errors="coerce")
row["%identity_min"] = vals.min()
row["%identity_max"] = vals.max()
if col_cov:
vals = pd.to_numeric(sub[col_cov], errors="coerce")
row["%coverage_min"] = vals.min()
row["%coverage_max"] = vals.max()
rows.append(row)
summary_iso_gene = pd.DataFrame(rows)
sum_iso_gene_tsv = "beta_lactam_summary_by_isolate_gene.tsv"
summary_iso_gene.to_csv(sum_iso_gene_tsv, sep='\t', index=False)
print(f">>> Isolate × gene summary written to: {sum_iso_gene_tsv}")
print(" (now includes 'Gene_symbol' and 'Sequence_name' annotation columns)")
# -------- optional Excel exports --------
try:
with pd.ExcelWriter("beta_lactam_all.xlsx") as xw:
df_beta.to_excel(xw, sheet_name="beta_lactam_all", index=False)
with pd.ExcelWriter("beta_lactam_summary.xlsx") as xw:
summary_gene.to_excel(xw, sheet_name="by_gene", index=False)
summary_iso_gene.to_excel(xw, sheet_name="by_isolate_gene", index=False)
print(">>> Excel workbooks written: beta_lactam_all.xlsx, beta_lactam_summary.xlsx")
except Exception as e:
print("WARNING: could not write Excel files:", e)
EOF
echo ">>> All done."
echo " - Individual reports: ${OUTDIR}/*.amrfinder.tsv"
echo " - Merged β-lactam table: beta_lactam_all.tsv"
echo " - Gene summary: beta_lactam_summary_by_gene.tsv (with isolate list)"
echo " - Isolate × gene summary: beta_lactam_summary_by_isolate_gene.tsv (with annotation)"
echo " - Excel (if pandas + openpyxl installed): beta_lactam_all.xlsx, beta_lactam_summary.xlsx"
3. Reporting and File Outputs
Files Generated:
beta_lactam_all.tsv: All β-lactam/beta-lactamase hits across genomes.beta_lactam_summary_by_gene.tsv: Per-gene summary, including a column with all isolate IDs.beta_lactam_summary_by_isolate_gene.tsv: Gene and isolate summary; includes “Gene symbol”, “Sequence name”, “Subclass”, annotation, and min/max identity/coverage.- If pandas is installed:
beta_lactam_all.xlsx, beta_lactam_summary.xlsx.
Description of improvements:
- Gene-level summary now lists isolates carrying each β-lactamase gene.
- Isolate × gene summary includes full annotation and quantitative metrics: gene symbol, sequence name, subclass, plus minimum and maximum percent identity/coverage.
Hand these files directly to collaborators:
beta_lactam_summary_by_isolate_gene.tsv or beta_lactam_summary.xlsx have all necessary gene and annotation information in final form.
4. Excel Export (if pandas is installed after the fact)
If the bacto environment lacks pandas, simply perform Excel conversion outside it:
mamba deactivate
python3 - << 'PYCODE'
import pandas as pd
df = pd.read_csv("beta_lactam_all.tsv", sep="\t")
df.to_excel("beta_lactam_all.xlsx", index=False)
print("Saved: beta_lactam_all.xlsx")
PYCODE
#Replace "," with ", " in beta_lactam_summary_by_gene.tsv so that the number can be correctly formatted; save it to a Excel-format.
mv beta_lactam_all.xlsx AMR_summary.xlsx # Then delete the first empty column.
mv beta_lactam_summary.xlsx BETA-LACTAM_summary.xlsx
mv beta_lactam_summary_by_gene.xlsx BETA-LACTAM_summary_by_gene.xlsx
Summary and Notes
- The system is fully automated from installation to reporting.
- All command lines are modular and suitable for direct inclusion in bioinformatics SOPs.
- Output files have expanded annotation and isolate information for downstream analytics and sharing.
- This approach ensures traceability, transparency, and rapid communication of β-lactamase annotation results for large datasets.
中文版本
基于NCBI AMRFinderPlus的自动化β-内酰胺酶注释流程
1. 安装与数据库设置
在bacto环境下安装AMRFinderPlus，并确保数据库更新：
mamba activate bacto
mamba install ncbi-amrfinderplus
mamba update ncbi-amrfinderplus
mamba activate bacto
amrfinder -u
- 这将自动下载最新数据库，并确保环境目录正确建立与软链接。
查询支持的物种选项：
amrfinder --list_organisms
2. 批量分析与脚本调用
使用如下脚本，高效批量筛查基因组β-内酰胺酶基因，并生成结果与汇总文件。
输入表格式：filename_TAB_organism，首行为表头。
cd ~/DATA/Data_Patricia_AMRFinderPlus_2025/genomes
./run_amrfinder_beta_lactam.sh genome_metadata.tsv
脚本逻辑简明，自动映射物种名、循环注释所有基因组、收集所有结果，之后调用Python脚本合并并筛选β-内酰胺酶基因。
3. 结果文件说明
beta_lactam_all.tsv：所有β-内酰胺酶相关基因注释全表beta_lactam_summary_by_gene.tsv：按基因注释，含所有分离株列表beta_lactam_summary_by_isolate_gene.tsv：分离株×基因详细表，包含注释、同源信息等- 若安装pandas：另有Excel版
beta_lactam_all.xlsx与beta_lactam_summary.xlsx
改进之处：
- 汇总表显式展示每个基因对应的分离株ID
- 分离株×基因表包含完整功能注释，identity/coverage等量化指标
直接将上述TSV或Excel表交给协作方即可，无需额外整理。
4. 补充Excel导出
如环境未装pandas，可以离线导出Excel：
mamba deactivate
python3 - << 'PYCODE'
import pandas as pd
df = pd.read_csv("beta_lactam_all.tsv", sep="\t")
df.to_excel("beta_lactam_all.xlsx", index=False)
print("Saved: beta_lactam_all.xlsx")
PYCODE
总结
- 步骤明确，可拓展与自动化。
- 输出表格式完善，满足批量汇报与协作需要。
- 所有命令与脚本可直接嵌入项目标准操作流程，支持可追溯和数据复用。
]
