宏基因组毒力因子分析：VFDB/VirulenceFinder¶

VFDB 毒力因子分类（14大类）：
├── Adherence（粘附）
├── Biofilm（生物膜）
├── Effector delivery system（效应物递送系统）
├── Exoenzyme（胞外酶）
├── Exotoxin（外毒素）
├── Immune evasion（免疫逃避）
├── Invasion（入侵）
├── Iron uptake（铁摄取）
├── Motility（运动性）
├── Nutritional/metabolic factor（营养/代谢因子）
├── Post-translational modification（翻译后修饰）
├── Regulation（调控）
├── Stress survival（应激存活）
└── Others（其他）

# ====== 方法1：使用DIAMOND比对VFDB ======

# 下载VFDB核心数据集（实验验证的VF）
wget http://www.mgc.ac.cn/VFs/Down/VFDB_setA_pro.fas.gz  # 下载核心蛋白序列
gunzip VFDB_setA_pro.fas.gz  # 解压

# 下载VFDB全集（包含预测的VF）
wget http://www.mgc.ac.cn/VFs/Down/VFDB_setB_pro.fas.gz  # 下载全集蛋白序列
gunzip VFDB_setB_pro.fas.gz

# 构建DIAMOND数据库
diamond makedb \
    --in VFDB_setA_pro.fas \  # 输入VFDB蛋白序列
    --db vfdb_diamond          # 输出数据库名

# 比对预测蛋白到VFDB
diamond blastp \
    --db vfdb_diamond \        # VFDB数据库
    --query predicted_proteins.faa \  # 输入预测的蛋白序列
    --out vfdb_results.tsv \   # 输出比对结果
    --outfmt 6 qseqid sseqid pident length evalue bitscore stitle \  # 输出格式
    --evalue 1e-10 \           # E-value阈值
    --id 80 \                  # 最低相似度80%
    --query-cover 70 \         # 最低查询覆盖度70%
    --max-target-seqs 1 \      # 每个查询只保留最佳匹配
    --threads 16               # 线程数

# ====== 方法2：使用ABRicate（支持多数据库）======

# 安装ABRicate
conda install -c bioconda abricate  # 安装abricate

# 查看可用数据库
abricate --list  # 列出所有可用数据库（包含vfdb）

# 运行ABRicate搜索VFDB
abricate \
    --db vfdb \                # 使用VFDB数据库
    contigs.fa \               # 输入组装的contigs
    --minid 80 \               # 最低相似度80%
    --mincov 70 \              # 最低覆盖度70%
    > abricate_vfdb_results.tsv  # 输出结果

# 汇总多个样本的结果
abricate --summary \
    sample1_results.tsv \      # 样本1结果
    sample2_results.tsv \      # 样本2结果
    > summary_matrix.tsv       # 输出汇总矩阵

# ====== 方法3：使用MetaVF工具包（2024年新工具）======
# MetaVF支持物种级别的VFG识别
# 安装和使用见 https://github.com/VFDB2/MetaVF

# 解析DIAMOND比对VFDB的结果
python3 << 'EOF'
import pandas as pd  # 导入pandas

# 读取比对结果
columns = ["qseqid", "sseqid", "pident", "length", "evalue", "bitscore", "stitle"]
df = pd.read_csv("vfdb_results.tsv", sep="\t", names=columns)  # 读取结果

# 从VFDB序列描述中解析信息
# VFDB格式示例：VFG000001(gb|AAA98077) (ace) [Adherence (VF0001)] [Vibrio cholerae]
import re
def parse_vfdb(stitle):
    """从VFDB描述行提取毒力因子信息"""
    vfg_match = re.search(r'(VFG\d+)', stitle)       # 提取VFG编号
    gene_match = re.search(r'\((\w+)\)', stitle)      # 提取基因名
    category_match = re.search(r'\[(.+?) \(', stitle)  # 提取分类
    species_match = re.findall(r'\[([^\]]+)\]', stitle)  # 提取物种

    return {
        "VFG_ID": vfg_match.group(1) if vfg_match else "NA",
        "gene": gene_match.group(1) if gene_match else "NA",
        "category": category_match.group(1) if category_match else "NA",
        "species": species_match[-1] if species_match else "NA"
    }

# 应用解析函数
parsed = df["stitle"].apply(parse_vfdb).apply(pd.Series)  # 解析每行
result = pd.concat([df, parsed], axis=1)  # 合并结果

# 统计各类别的VF数量
print("=== 毒力因子类别分布 ===")
print(result["category"].value_counts())  # 打印各类别数量

# 统计各物种的VF数量
print("\n=== 携带VF的物种分布 ===")
print(result["species"].value_counts().head(10))  # 打印Top10物种
EOF

# 在线版本（推荐初学者使用）：
# https://cge.food.dtu.dk/services/VirulenceFinder/
# 支持上传FASTA文件，快速获得结果

# 本地安装
conda install -c bioconda virulencefinder  # 安装virulencefinder

# 下载数据库
git clone https://bitbucket.org/genomicepidemiology/virulencefinder_db.git

# 运行VirulenceFinder
python -m virulencefinder \
    -i contigs.fa \                    # 输入序列文件
    -o virulencefinder_output/ \       # 输出目录
    -p virulencefinder_db/ \           # 数据库路径
    -s "Escherichia coli" \            # 指定物种
    --min_cov 0.6 \                    # 最低覆盖度60%
    --threshold 0.9                    # 最低相似度90%

# PathoFact同时预测VF和ARG
# 安装（需要Snakemake）
git clone https://git-r3lab.uni.lu/laura.zinnliebenstein/pathoFact.git  # 克隆仓库

# 运行PathoFact
snakemake \
    --snakefile PathoFact/Snakefile \  # 指定Snakefile
    --configfile config.yaml \         # 配置文件
    --cores 16                         # 核心数

# 在你的T2D研究中分析VF的示例代码
# 第一步：比对MAG蛋白到VFDB
diamond blastp --db vfdb_diamond --query t2d_mag_proteins.faa \
    --out t2d_vf_results.tsv --outfmt 6 --evalue 1e-10 --id 80 --threads 16

# 第二步：比较T2D组和对照组的VF差异
python3 << 'EOF'
import pandas as pd
from scipy import stats  # 导入统计模块

# 假设已经有每个样本的VF丰度矩阵
# vf_abundance: 行=样本, 列=VF基因, 值=丰度
# metadata: 包含Group列(T2D/Control)

# 对每个VF基因做Mann-Whitney U检验
for vf_gene in vf_abundance.columns:
    t2d_values = vf_abundance.loc[metadata["Group"]=="T2D", vf_gene]  # T2D组
    ctrl_values = vf_abundance.loc[metadata["Group"]=="Control", vf_gene]  # 对照组
    stat, pval = stats.mannwhitneyu(t2d_values, ctrl_values)  # U检验
    if pval < 0.05:
        print(f"{vf_gene}: p={pval:.4f}")  # 输出显著差异的VF
EOF

# 工具选择
最全面的VF分析 → VFDB + DIAMOND比对
快速多数据库筛查 → ABRicate（支持VFDB/VirulenceFinder等）
在线快速检测 → VirulenceFinder网页版
宏基因组物种级分析 → MetaVF（2024年新工具）
VF + ARG联合分析 → PathoFact

# VFDB数据集选择
setA = 核心数据集（实验验证的VF，推荐）
setB = 全集（包含预测的VF同源序列）

# 比对阈值推荐
标准分析：identity ≥ 80%, coverage ≥ 70%, e-value ≤ 1e-10
严格分析：identity ≥ 90%, coverage ≥ 80%, e-value ≤ 1e-20

# 14大类毒力因子速记
粘附、生物膜、分泌系统、胞外酶、外毒素、免疫逃避、入侵
铁摄取、运动性、营养因子、翻译后修饰、调控、应激存活、其他