宏基因组抗性基因分型¶

# CARD数据库结构
# ARO (Antibiotic Resistance Ontology) - 基因分类系统
# - Drug class: 抗药物类别
# - Resistance mechanism: 耐药机制
# - AMR gene family: 基因家族
# - Target: 药物靶点

# 下载最新CARD数据库（当前数据版本4.0.1，2025年2月更新）
wget https://card.mcmaster.ca/latest/data
tar -xvf data ./card.json

# 安装RGI工具（当前版本6.0.5，推荐conda安装）
# conda install -c bioconda rgi  # 推荐方式
pip install rgi  # 或pip安装
rgi auto_load --card_json card.json

# === 宏基因组预处理 ===

# 1. 质控
fastp -i sample_R1.fq.gz -I sample_R2.fq.gz \
    -o clean_R1.fq.gz -O clean_R2.fq.gz \
    --detect_adapter_for_pe --thread 16

# 2. 去宿主序列（人类样本）
bowtie2 -x human_genome -1 clean_R1.fq.gz -2 clean_R2.fq.gz \
    --threads 16 --very-sensitive \
    --un-conc-gz nonhost_R%.fq.gz \
    -S /dev/null

# 3. 宏基因组组装（可选，用于contig-level分析）
megahit -1 nonhost_R1.fq.gz -2 nonhost_R2.fq.gz \
    -o megahit_assembly --min-contig-len 500 \
    --num-cpu-threads 32 --memory 0.9

# 4. 基因预测（组装后）
prodigal -i megahit_assembly/final.contigs.fa \
    -o genes.gff -a proteins.faa -d genes.fna \
    -p meta -f gff

# === RGI 分析 ===

# 1. 基于Contig的RGI分析（核酸序列）
rgi main \
    --input_sequence megahit_assembly/final.contigs.fa \
    --output_file sample_rgi \
    --input_type contig \
    --alignment_tool BLAST \
    --num_threads 16 \
    --clean \
    --include_loose

# 2. 基于蛋白序列的RGI分析
rgi main \
    --input_sequence proteins.faa \
    --output_file sample_rgi_protein \
    --input_type protein \
    --alignment_tool DIAMOND \
    --num_threads 16

# 3. 基于Reads的RGI bwt分析（不需组装）
rgi bwt \
    --read_one nonhost_R1.fq.gz \
    --read_two nonhost_R2.fq.gz \
    --output_file sample_rgi_bwt \
    --aligner bowtie2 \
    --threads 16 \
    --include_wildcard

# 4. 结果解析
# 输出文件: sample_rgi.txt (tab-delimited)
# 关键列: Best_Hit_ARO, Drug Class, Resistance Mechanism,
#         AMR Gene Family, Cut_Off (Perfect/Strict/Loose)

# RGI结果解析与统计
import pandas as pd
import matplotlib.pyplot as plt

# 读取RGI结果
rgi_results = pd.read_csv("sample_rgi.txt", sep="\t")

# 过滤只保留 Perfect + Strict hits
confident_hits = rgi_results[rgi_results['Cut_Off'].isin(['Perfect', 'Strict'])]

# 按Drug Class统计
drug_class_counts = confident_hits['Drug Class'].str.split('; ').explode().value_counts()
print("Top 10 Drug Classes:")
print(drug_class_counts.head(10))

# 按Resistance Mechanism统计
mech_counts = confident_hits['Resistance Mechanism'].value_counts()
print("\nResistance Mechanisms:")
print(mech_counts)

# 按AMR Gene Family统计
family_counts = confident_hits['AMR Gene Family'].value_counts()
print("\nTop 10 AMR Gene Families:")
print(family_counts.head(10))

# === StarAMR 安装与运行 ===

# 安装
pip install staramr
staramr db build --dir staramr_db

# 运行（输入为组装的contigs或完整基因组）
staramr search \
    --database staramr_db \
    --output-dir staramr_results \
    --pid-threshold 90 \
    --percent-length-overlap-resfinder 60 \
    --nprocs 16 \
    assembly1.fasta assembly2.fasta assembly3.fasta

# 输出文件说明
# resfinder.tsv     - 获得性ARG检测结果
# pointfinder.tsv   - 点突变耐药
# plasmidfinder.tsv - 质粒replicon
# mlst.tsv          - MLST分型
# summary.tsv       - 综合汇总
# detailed_summary.tsv - 详细汇总

# 查看结果
cat staramr_results/summary.tsv

# === ResFinder 命令行使用（当前版本4.7.2，2025年4月更新）===
pip install resfinder

# 运行 ResFinder
python -m resfinder \
    --inputfasta assembly.fasta \
    --outputPath resfinder_output \
    --acquired \
    --species "Escherichia coli" \
    --min_cov 0.6 \
    --threshold 0.9 \
    --db_path_res resfinder_db

# 同时检测点突变（需指定物种）
python -m resfinder \
    --inputfasta assembly.fasta \
    --outputPath resfinder_output \
    --acquired --point \
    --species "Escherichia coli" \
    --db_path_res resfinder_db \
    --db_path_point pointfinder_db

# ResFinder 也支持reads输入
python -m resfinder \
    --inputfastq reads_R1.fq.gz reads_R2.fq.gz \
    --outputPath resfinder_reads_output \
    --acquired \
    --species "Klebsiella pneumoniae"

# === ARG定量（reads-based）===

# 方法1：比对到CARD序列数据库
# 构建比对索引
bowtie2-build card_nucleotide_seqs.fasta card_index

# 比对
bowtie2 -x card_index \
    -1 nonhost_R1.fq.gz -2 nonhost_R2.fq.gz \
    --threads 16 --very-sensitive \
    -S sample_vs_card.sam

# 计算覆盖度和丰度
samtools view -bS -F 4 sample_vs_card.sam | samtools sort > sample_vs_card.sorted.bam
samtools index sample_vs_card.sorted.bam
samtools idxstats sample_vs_card.sorted.bam > sample_card_counts.txt

# 方法2：使用ARGs-OAP流程（推荐宏基因组）
# ARGs-OAP = Antibiotic Resistance Genes Online Analysis Pipeline
# https://github.com/xiaole99/ARGs-OAP-v3.0

# 方法3：ShortBRED（高特异性短标记序列）
shortbred_identify.py \
    --goi card_proteins.faa \
    --ref uniref90.fasta \
    --markers card_markers.faa \
    --tmp tmp_identify

shortbred_quantify.py \
    --markers card_markers.faa \
    --wgs nonhost_R1.fq.gz \
    --results sample_arg_shortbred.tsv \
    --tmp tmp_quantify \
    --threads 16

# ARG丰度归一化
import pandas as pd
import numpy as np

# 读取比对计数
counts = pd.read_csv("sample_card_counts.txt", sep="\t",
                      header=None, names=['gene', 'length', 'mapped', 'unmapped'])
counts = counts[counts['gene'] != '*']

# 总reads数
total_reads = counts['mapped'].sum() + counts['unmapped'].sum()

# 计算RPKM
counts['RPKM'] = (counts['mapped'] * 1e9) / (counts['length'] * total_reads)

# 计算copies per 16S
# 假设已知16S比对reads数
reads_16s = 50000  # 比对到16S的reads数
gene_16s_length = 1500  # 16S长度

copies_16s = (reads_16s * 1e3) / gene_16s_length
counts['copies_per_16S'] = (counts['mapped'] / counts['length'] * 1e3) / copies_16s

# 按药物类别汇总
# 需要ARO注释文件关联gene -> drug class
aro_mapping = pd.read_csv("aro_categories.tsv", sep="\t")
merged = counts.merge(aro_mapping, on='gene')
drug_class_abundance = merged.groupby('Drug Class')['RPKM'].sum().sort_values(ascending=False)

# === R语言ARG可视化 ===
library(ggplot2)
library(pheatmap)
library(vegan)

# 读取多样本ARG丰度矩阵
arg_matrix <- read.csv("arg_abundance_matrix.csv", row.names = 1)
# 行=ARG, 列=样本

# 热图展示ARG丰度
pheatmap(log2(arg_matrix + 1),
         clustering_distance_rows = "euclidean",
         clustering_distance_cols = "euclidean",
         show_rownames = FALSE,
         annotation_col = sample_meta)

# ARG多样性分析（与微生物多样性类似）
arg_diversity <- diversity(t(arg_matrix), index = "shannon")
arg_richness <- specnumber(t(arg_matrix))

# PCoA分析
arg_dist <- vegdist(t(arg_matrix), method = "bray")
pcoa_result <- cmdscale(arg_dist, k = 2, eig = TRUE)
pcoa_df <- data.frame(
  PC1 = pcoa_result$points[, 1],
  PC2 = pcoa_result$points[, 2],
  Group = sample_meta$group
)
ggplot(pcoa_df, aes(x = PC1, y = PC2, color = Group)) +
  geom_point(size = 3) +
  stat_ellipse() +
  theme_minimal()

# PERMANOVA检验组间差异
adonis2(arg_dist ~ Group, data = sample_meta)

# === 快速ARG分型流水线 ===

# 单基因组/分离株
rgi main -i genome.fasta -o rgi_result -t contig -a DIAMOND -n 8
staramr search -o star_result genome.fasta

# 宏基因组（基于reads）
rgi bwt --read_one R1.fq.gz --read_two R2.fq.gz \
    --output_file meta_rgi --aligner bowtie2 --threads 16

# 宏基因组（基于组装）
megahit -1 R1.fq.gz -2 R2.fq.gz -o assembly
prodigal -i assembly/final.contigs.fa -a proteins.faa -p meta
rgi main -i proteins.faa -o assembly_rgi -t protein -a DIAMOND -n 16

# 数据库更新
rgi auto_load --card_json /path/to/card.json
staramr db build --dir /path/to/staramr_db