669 细菌耐药机制与基因组¶

# AMRFinderPlus安装
conda install -c bioconda ncbi-amrfinderplus  # 安装
amrfinder --update  # 更新数据库

# 从基因组检测耐药基因
amrfinder \
  -n genome.fasta \           # 核酸序列输入
  -p proteins.faa \           # 蛋白质序列输入（可选，提高准确性）
  -g gff_annotation.gff \     # GFF注释文件
  -O Escherichia \            # 物种（用于点突变检测）
  --plus \                    # 额外检测毒力因子和应激基因
  -o amrfinder_results.tsv \  # 输出结果
  --threads 8                 # 线程数

# 输出列包含：
# Gene symbol, Sequence name, Element type, Element subtype,
# Class, Subclass, Method, Target length, % Coverage, % Identity

# RGI安装
conda install -c bioconda rgi  # 安装
rgi load --card_json card.json  # 加载CARD数据库

# 从基因组检测耐药基因
rgi main \
  --input_sequence genome.fasta \  # 输入基因组
  --output_file rgi_results \      # 输出前缀
  --input_type contig \            # 输入类型：contig
  --alignment_tool DIAMOND \       # 比对工具
  --num_threads 8 \                # 线程数
  --clean                          # 清理临时文件

# 从宏基因组reads检测
rgi bwt \
  --read_one reads_R1.fastq.gz \   # 正向reads
  --read_two reads_R2.fastq.gz \   # 反向reads
  --output_file rgi_meta \         # 输出前缀
  --aligner bowtie2 \              # 比对器
  --threads 8                      # 线程数

# RGI检测模型类型：
# Protein Homolog Model: 基于蛋白同源性
# Protein Variant Model: 基于已知耐药突变
# rRNA Variant Model: rRNA突变检测
# Protein Overexpression Model: 过表达检测

# ResFinder安装
conda install -c bioconda resfinder  # 安装

# 检测获得性耐药基因
run_resfinder.py \
  -ifq reads_R1.fastq.gz reads_R2.fastq.gz \  # 输入reads
  -o resfinder_results/ \    # 输出目录
  -s "Escherichia coli" \    # 物种（用于点突变）
  --acquired \               # 检测获得性耐药基因
  --point \                  # 检测点突变
  -db_res /path/to/resfinder_db \  # 耐药基因数据库
  -db_point /path/to/pointfinder_db  # 点突变数据库

# 解析和比较AMR检测结果
import pandas as pd  # 数据处理
import matplotlib.pyplot as plt  # 绑图
import seaborn as sns  # 热图

# 1. 读取AMRFinderPlus结果
amr = pd.read_csv("amrfinder_results.tsv", sep='\t')
print(f"检测到 {len(amr)} 个耐药/毒力基因")

# 按抗生素类别统计
class_counts = amr["Class"].value_counts()
print("\n抗生素类别分布:")
print(class_counts)

# 2. 构建样本-耐药基因矩阵（多样本比较）
import os
import glob

all_results = []
for f in glob.glob("amr_results/*.tsv"):
    sample = os.path.basename(f).replace("_amr.tsv", "")
    df = pd.read_csv(f, sep='\t')
    df["sample"] = sample
    all_results.append(df)

combined = pd.concat(all_results)

# 构建存在/缺失矩阵
amr_matrix = combined.pivot_table(
    index="sample",
    columns="Gene symbol",
    values="% Identity of matching region",
    aggfunc="max"
).fillna(0)
amr_binary = (amr_matrix > 0).astype(int)  # 二值化

# 3. 热图可视化
fig, ax = plt.subplots(figsize=(16, 8))
sns.heatmap(amr_binary,
            cmap="YlOrRd",  # 颜色方案
            linewidths=0.5,
            cbar_kws={"label": "存在(1)/缺失(0)"},
            ax=ax)
ax.set_title("耐药基因存在/缺失热图")
ax.set_ylabel("样本")
ax.set_xlabel("耐药基因")
plt.tight_layout()
plt.savefig("amr_heatmap.png", dpi=150)

# 4. 基因型-表型比较
phenotype = pd.read_csv("mic_results.csv", index_col=0)  # MIC检测结果
# 计算基因型预测的敏感性和特异性
def genotype_phenotype_concordance(amr_binary, phenotype, drug_gene_map):
    """计算基因型预测表型的准确性"""
    for drug, genes in drug_gene_map.items():
        present_genes = [g for g in genes if g in amr_binary.columns]
        if not present_genes or drug not in phenotype.columns:
            continue
        # 基因型预测：任意一个相关基因存在→预测耐药
        geno_pred = amr_binary[present_genes].max(axis=1)
        pheno_true = (phenotype[drug] == "R").astype(int)

        # 计算准确性
        common = geno_pred.index.intersection(pheno_true.index)
        tp = ((geno_pred[common] == 1) & (pheno_true[common] == 1)).sum()
        tn = ((geno_pred[common] == 0) & (pheno_true[common] == 0)).sum()
        fp = ((geno_pred[common] == 1) & (pheno_true[common] == 0)).sum()
        fn = ((geno_pred[common] == 0) & (pheno_true[common] == 1)).sum()

        sens = tp / (tp + fn) if (tp + fn) > 0 else 0
        spec = tn / (tn + fp) if (tn + fp) > 0 else 0
        print(f"{drug}: 灵敏度={sens:.2f}, 特异度={spec:.2f}")

# 从宏基因组数据检测耐药基因
# 方法1: ShortBRED（marker-based）
shortbred_quantify.py \
  --markers ShortBRED_CARD_2024_markers.faa \  # CARD标记基因
  --wgs metagenome.fastq.gz \  # 宏基因组reads
  --results shortbred_results.tsv \  # 输出
  --threads 8

# 方法2: AMR++管线
# nextflow run https://github.com/Microbial-Ecology-Group/AMRplusplus
nextflow run amrplusplus \
  --reads "reads/*_{1,2}.fastq.gz" \  # 输入reads
  --output amrpp_results/ \           # 输出
  --pipeline standard                 # 标准管线

# AMR基因组分析流程
基因组/MAG/reads → 耐药基因检测
  AMRFinderPlus: NCBI开发，含毒力基因
  RGI (CARD): 最全面，含突变模型
  ResFinder: 获得性基因+点突变
  ABRicate: 快速筛查，支持多数据库
  → 构建AMR profile → 基因型-表型对比

# 数据库对比
CARD: 8582本体术语，6442参考序列，最全面
NDARO/AMRFinderPlus: 实时更新，侧重获得性基因
ResFinder: Enterobacteriaceae优化
MEGARes: 宏基因组优化

# 工具选择
精确分析: AMRFinderPlus + RGI (推荐组合)
快速筛查: ABRicate
宏基因组: RGI bwt / AMR++ / ShortBRED
多物种覆盖(2025): AmrProfiler (18,608物种)

# 临床报告关键信息
1. 检测到的耐药基因列表
2. 对应的抗生素耐药预测
3. 基因位置(染色体/质粒)
4. 可移动性评估(是否可水平转移)