宏基因组抗生素耐药基因分析：CARD/ResFinder/AMRFinderPlus¶

# 安装RGI（CARD的分析工具）
conda install -c bioconda rgi  # 用conda安装rgi

# 下载最新CARD数据库
wget https://card.mcmaster.ca/latest/data  # 下载数据库
tar -xvf data ./card.json                  # 解压

# 加载数据库
rgi load \
    --card_json card.json \      # 加载CARD数据库
    --local                      # 保存到本地

# ====== 模式1：分析基因组/MAG（蛋白序列）======
rgi main \
    --input_sequence predicted_proteins.faa \  # 输入蛋白序列
    --output_file rgi_genome_results \         # 输出前缀
    --input_type protein \                     # 输入类型：蛋白质
    --local \                                  # 使用本地数据库
    --clean \                                  # 清理临时文件
    -n 16                                      # 线程数

# ====== 模式2：分析宏基因组reads ======
rgi bwt \
    --read_one sample_R1.fastq.gz \  # 输入正向reads
    --read_two sample_R2.fastq.gz \  # 输入反向reads
    --output_file rgi_meta_results \ # 输出前缀
    --local \                        # 使用本地数据库
    --aligner bowtie2 \              # 比对工具：bowtie2或kma
    -n 16                            # 线程数

# ====== 模式3：分析组装后的contigs ======
rgi main \
    --input_sequence contigs.fa \    # 输入组装的contigs
    --output_file rgi_contig_results \  # 输出前缀
    --input_type contig \            # 输入类型：contig
    --local \
    --clean \
    -n 16

# 查看结果
cat rgi_genome_results.txt  # 查看结果表格

# 关键列说明：
# ORF_ID          — 开放阅读框ID
# Best_Hit_ARO    — 最佳匹配的ARO术语（耐药基因名称）
# Drug Class      — 耐药的药物类别
# Resistance Mechanism — 耐药机制
# AMR Gene Family — ARG基因家族
# % Identity      — 序列相似度百分比
# Model_type      — 检测模型类型：
#   Perfect: 与已知ARG 100%匹配
#   Strict: 高置信度匹配（推荐使用的阈值）
#   Loose: 低置信度匹配（可能是假阳性）

# 安装AMRFinderPlus
conda install -c bioconda ncbi-amrfinderplus  # 安装

# 更新数据库（每次使用前建议更新）
amrfinder --update  # 更新到最新数据库版本

# ====== 分析蛋白序列 ======
amrfinder \
    --protein predicted_proteins.faa \  # 输入蛋白序列
    --output amrfinder_results.tsv \   # 输出结果文件
    --threads 16 \                     # 线程数
    --plus                             # 开启plus模式（包含应激/毒力基因）

# ====== 分析核酸序列 ======
amrfinder \
    --nucleotide contigs.fa \          # 输入核酸序列
    --output amrfinder_nuc_results.tsv \  # 输出文件
    --threads 16 \
    --plus

# ====== 同时分析蛋白+核酸+GFF注释（最准确）======
amrfinder \
    --protein predicted_proteins.faa \  # 蛋白序列
    --nucleotide contigs.fa \           # 核酸序列
    --gff predicted_genes.gff \         # GFF注释文件
    --output amrfinder_full.tsv \       # 输出文件
    --organism Escherichia \            # 指定物种（用于检测点突变）
    --threads 16 \
    --plus

# 关键列说明：
# Gene symbol     — 基因名称
# Sequence name   — 所在序列ID
# Element type    — AMR/STRESS/VIRULENCE
# Element subtype — POINT（点突变）/AMR（获得性基因）
# Class           — 抗生素类别
# Subclass        — 抗生素亚类
# Method          — 检测方法：EXACTX/ALLELEX/BLASTX/HMM/POINTX
# % Coverage      — 覆盖度
# % Identity      — 相似度

# 安装ResFinder本地版
conda install -c bioconda resfinder  # 安装

# 下载数据库
git clone https://bitbucket.org/genomicepidemiology/resfinder_db.git  # 下载耐药基因数据库
git clone https://bitbucket.org/genomicepidemiology/pointfinder_db.git  # 下载点突变数据库

# 运行ResFinder
python -m resfinder \
    -ifq sample_R1.fastq.gz sample_R2.fastq.gz \  # 输入reads
    -o resfinder_output/ \   # 输出目录
    -db_res resfinder_db/ \  # 耐药基因数据库路径
    -db_point pointfinder_db/ \  # 点突变数据库路径
    -s "Escherichia coli" \  # 指定物种（用于点突变检测）
    --min_cov 0.6 \          # 最低覆盖度60%
    --threshold 0.9          # 最低相似度90%

# 也可以使用在线版本：
# https://cge.food.dtu.dk/services/ResFinder/

# ====== 推荐的完整分析流程 ======

# 第一步：质控和去宿主
fastp -i R1.fq.gz -I R2.fq.gz -o clean_R1.fq.gz -O clean_R2.fq.gz  # 质控
bowtie2 -x human_ref -1 clean_R1.fq.gz -2 clean_R2.fq.gz \
    --un-conc-gz non_host_%.fq.gz -S /dev/null  # 去宿主reads

# 第二步：组装（可选，提高灵敏度）
megahit -1 non_host_1.fq.gz -2 non_host_2.fq.gz -o assembly/ -t 16  # 组装

# 第三步：基因预测
prodigal -i assembly/final.contigs.fa \
    -o genes.gff \           # 输出GFF注释
    -a proteins.faa \        # 输出蛋白序列
    -p meta                  # 宏基因组模式

# 第四步：多工具并行分析ARG
# 4a. CARD/RGI分析
rgi main -i proteins.faa -o rgi_results -t protein --local -n 16

# 4b. AMRFinderPlus分析
amrfinder --protein proteins.faa --nucleotide assembly/final.contigs.fa \
    --gff genes.gff -o amrfinder_results.tsv -t 16 --plus

# 第五步：结果整合与统计
python3 << 'EOF'
import pandas as pd  # 导入pandas

# 读取RGI结果
rgi = pd.read_csv("rgi_results.txt", sep="\t")  # 读取RGI结果
rgi_strict = rgi[rgi["Cut_Off"] == "Strict"]     # 只保留Strict级别的结果

# 读取AMRFinderPlus结果
amr = pd.read_csv("amrfinder_results.tsv", sep="\t")  # 读取AMRFinder结果
amr_genes = amr[amr["Element type"] == "AMR"]          # 只保留AMR类型

# 统计各药物类别的ARG数量
print("=== CARD/RGI 结果 ===")
print(rgi_strict["Drug Class"].value_counts().head(10))  # 打印Top10药物类别

print("\n=== AMRFinderPlus 结果 ===")
print(amr_genes["Class"].value_counts().head(10))  # 打印Top10药物类别

# 统计耐药机制分布
print("\n=== 耐药机制分布 ===")
print(rgi_strict["Resistance Mechanism"].value_counts())  # 各机制的ARG数量
EOF

# 工具选择
研究级全面分析 → CARD/RGI（最全面的数据库）
临床提交/公卫监测 → AMRFinderPlus（NCBI官方）
快速在线筛查 → ResFinder（DTU在线平台）
探索新型ARG → ResFinderFG v2.0（功能宏基因组来源）
深度学习检测新ARG → DeepARG

# CARD/RGI结果筛选标准
Perfect: 100%匹配已知ARG → 高度确信
Strict: 高置信度比特分数阈值以上 → 推荐使用
Loose: 低置信度 → 通常不用，容易假阳性

# 常见耐药机制
antibiotic inactivation — 抗生素灭活（如β-内酰胺酶）
antibiotic target alteration — 靶标修饰
antibiotic efflux — 外排泵
reduced permeability — 降低通透性
antibiotic target replacement — 靶标替换
antibiotic target protection — 靶标保护

# 关键命令速查
rgi main -i proteins.faa -o out -t protein --local -n 16
amrfinder --protein prot.faa -o out.tsv -t 16 --plus
amrfinder --update