734. 宏基因组酶活性预测¶
一句话概述：从宏基因组数据中预测微生物群落能产生哪些酶、酶活性有多强——就像拿到一个工厂的设备清单，推算这个工厂能生产什么、产能多大。
核心知识点速查表¶
| 概念 | 白话解释 | 关键工具 |
|---|---|---|
| EC号 | 酶的"身份证号"，4级分类（如EC 3.2.1.1=淀粉酶） | KEGG, BRENDA |
| CAZyme | 碳水化合物活性酶，降解/修饰糖类的酶家族 | dbCAN, CAZy |
| 功能丰度 | 某个酶在群落中有多少拷贝（基因数） | HUMAnN3 |
| 酶活预测 | 从序列推测蛋白是否真的有酶活性 | DeepEC, CLEAN |
| 底物预测 | 预测酶催化的具体底物是什么 | ESP, EnzymeMiner |
| 最适温度/pH | 预测酶的最佳工作条件 | Tm预测模型 |
一、背景与原理（白话版）¶
1.1 为什么需要预测酶活性？¶
宏基因组告诉你群落里有什么基因，但： - 有基因 ≠ 有蛋白（可能不表达） - 有蛋白 ≠ 有活性（可能是假基因） - 有活性 ≠ 活性强（可能效率很低）
酶活性预测就是在没有实验条件下，用计算方法估计： 1. 这个基因编码的蛋白是不是真的酶 2. 如果是，它的催化效率如何 3. 它的底物特异性是什么
1.2 应用场景¶
| 场景 | 具体问题 | 分析方法 |
|---|---|---|
| 肠道菌群 | 哪些菌能降解膳食纤维？ | CAZyme分析 |
| 工业酶发现 | 宏基因组中有没有新的耐热蛋白酶？ | 功能筛选+Tm预测 |
| 抗生素耐药 | 哪些菌携带β-内酰胺酶？ | 耐药基因+酶活预测 |
| 环境修复 | 降解污染物的酶在哪些菌里？ | 功能基因检索 |
二、CAZyme分析流程（最常用）¶
2.1 dbCAN3安装与使用¶
# ===== 安装dbCAN3（CAZyme注释的金标准工具）=====
conda create -n cazyme python=3.9 # 创建专用环境
conda activate cazyme # 激活环境
# 安装dbCAN
pip install dbcan # 安装dbCAN主程序
# 下载数据库
dbcan_build --cpus 8 --db-dir db # 自动下载并构建数据库（约10GB）
# 如果自动下载失败，手动下载
mkdir -p db && cd db # 创建并进入数据库目录
wget https://bcb.unl.edu/dbCAN2/download/dbCAN-HMMdb-V12.txt # HMM数据库
wget https://bcb.unl.edu/dbCAN2/download/CAZyDB.fa # 序列数据库
wget https://bcb.unl.edu/dbCAN2/download/tcdb.fa # 转运蛋白数据库
wget https://bcb.unl.edu/dbCAN2/download/stp.hmm # 信号肽数据库
hmmpress dbCAN-HMMdb-V12.txt # 构建HMM索引
diamond makedb --in CAZyDB.fa -d CAZy # 构建DIAMOND索引
2.2 从宏基因组预测CAZyme¶
# ===== 第一步：基因预测 =====
# 从宏基因组组装结果中预测蛋白编码基因
prodigal -i assembly.fa \ # 输入：宏基因组组装结果
-a proteins.faa \ # 输出：蛋白序列
-d genes.fna \ # 输出：核酸序列
-o genes.gff \ # 输出：GFF注释
-p meta \ # 宏基因组模式
-f gff # 输出格式为GFF
# ===== 第二步：运行dbCAN注释 =====
run_dbcan proteins.faa protein \ # 输入蛋白序列，指定类型为protein
--db_dir db/ \ # 数据库目录
--out_dir dbcan_output/ \ # 输出目录
--dia_cpu 8 \ # DIAMOND使用的CPU数
--hmm_cpu 8 \ # HMMER使用的CPU数
--tools all # 使用所有三种方法（HMMER+DIAMOND+eCAMI）
# ===== 结果文件说明 =====
ls dbcan_output/
# overview.txt # 汇总结果（推荐用这个）
# diamond.out # DIAMOND比对结果
# hmmer.out # HMMER搜索结果
# eCAMI.out # eCAMI深度学习结果
# substrate.out # 底物预测结果（新功能）
2.3 结果解析¶
# ===== 解析dbCAN结果 =====
import pandas as pd # 导入pandas
# 读取overview结果
df = pd.read_csv("dbcan_output/overview.txt", sep="\t") # 读取汇总表
print(df.head()) # 查看前几行
# Gene_ID HMMER eCAMI DIAMOND #ofTools
# gene_1 GH13 GH13 GH13 3 ← 三种方法都支持，非常可靠
# gene_2 GH5 - GH5 2 ← 两种方法支持，较可靠
# gene_3 - GT2 - 1 ← 只有一种方法支持，不太可靠
# 过滤：至少2种方法一致才保留
reliable = df[df["#ofTools"] >= 2] # 保留至少两种方法支持的结果
print(f"可靠的CAZyme基因数: {len(reliable)}") # 打印数量
# 统计CAZyme家族分布
# 提取CAZyme家族信息（取HMMER列作为主注释）
cazyme_counts = reliable["HMMER"].value_counts() # 统计每个家族的基因数
print("Top 20 CAZyme families:")
print(cazyme_counts.head(20)) # 打印前20个最多的家族
# ===== 按CAZyme大类统计 =====
def get_class(family_str):
"""从CAZyme家族名提取大类"""
if pd.isna(family_str) or family_str == "-": # 处理缺失值
return "Unknown"
family = family_str.split("+")[0].split("_")[0] # 取第一个家族
if family.startswith("GH"): # 糖苷水解酶
return "GH (Glycoside Hydrolases)"
elif family.startswith("GT"): # 糖基转移酶
return "GT (Glycosyl Transferases)"
elif family.startswith("PL"): # 多糖裂解酶
return "PL (Polysaccharide Lyases)"
elif family.startswith("CE"): # 碳水化合物酯酶
return "CE (Carbohydrate Esterases)"
elif family.startswith("AA"): # 辅助活性酶
return "AA (Auxiliary Activities)"
elif family.startswith("CBM"): # 碳水化合物结合模块
return "CBM (Carbohydrate-Binding Modules)"
else:
return "Other"
reliable["Class"] = reliable["HMMER"].apply(get_class) # 添加大类列
class_counts = reliable["Class"].value_counts() # 按大类统计
# 绘制饼图
import matplotlib.pyplot as plt # 导入绑图库
plt.figure(figsize=(10, 8)) # 设置画布
colors = ['#e74c3c', '#3498db', '#2ecc71', '#f39c12', '#9b59b6', '#1abc9c', '#95a5a6']
plt.pie(class_counts.values, labels=class_counts.index, # 绘制饼图
colors=colors, autopct='%1.1f%%', startangle=140)
plt.title("CAZyme Class Distribution in Metagenome") # 标题
plt.tight_layout()
plt.savefig("cazyme_distribution.png", dpi=300) # 保存
plt.show()
三、DeepEC：深度学习酶分类¶
# ===== 安装DeepEC =====
# DeepEC用深度学习从蛋白序列预测EC号
git clone https://github.com/deepec/DeepEC.git # 克隆仓库
cd DeepEC # 进入目录
pip install -r requirements.txt # 安装依赖
# 运行DeepEC预测
python deepec.py \
-i proteins.faa \ # 输入：蛋白序列文件
-o deepec_output/ \ # 输出目录
--device cpu # 使用CPU（有GPU可改为cuda:0）
# 结果文件
# deepec_output/
# ├── DeepEC_Result.txt # EC号预测结果
# └── DeepEC_Result.log # 运行日志
# ===== 解析DeepEC结果 =====
import pandas as pd # 导入pandas
# 读取结果
results = pd.read_csv("deepec_output/DeepEC_Result.txt", sep="\t",
names=["Protein_ID", "EC_Number", "Score"]) # 自定义列名
# 过滤高置信度预测
high_conf = results[results["Score"] >= 0.9] # 只保留置信度>=0.9的预测
print(f"高置信度酶预测数: {len(high_conf)}")
# 统计EC大类分布
# EC号格式：EC X.X.X.X
# 第一级：1=氧化还原酶 2=转移酶 3=水解酶 4=裂合酶 5=异构酶 6=连接酶 7=转位酶
high_conf["EC_Class"] = high_conf["EC_Number"].str.split(".").str[0] # 提取第一级分类
ec_class_names = {
"1": "Oxidoreductases (氧化还原酶)",
"2": "Transferases (转移酶)",
"3": "Hydrolases (水解酶)",
"4": "Lyases (裂合酶)",
"5": "Isomerases (异构酶)",
"6": "Ligases (连接酶)",
"7": "Translocases (转位酶)"
}
high_conf["EC_Class_Name"] = high_conf["EC_Class"].map(ec_class_names) # 映射中文名
print(high_conf["EC_Class_Name"].value_counts()) # 打印各类酶的数量
四、HUMAnN3功能丰度定量¶
# ===== 用HUMAnN3定量酶的丰度 =====
# HUMAnN3可以从宏基因组reads直接定量代谢通路和酶的丰度
# 安装
conda install -c bioconda humann # 安装HUMAnN3
# 运行
humann --input sample_R1.fastq.gz \ # 输入：质控后的reads
--output humann_output/ \ # 输出目录
--threads 8 \ # 线程数
--search-mode uniref90 # 搜索模式（uniref90精度高）
# 结果文件
# humann_output/
# ├── sample_genefamilies.tsv # 基因家族丰度（UniRef90）
# ├── sample_pathabundance.tsv # 代谢通路丰度
# └── sample_pathcoverage.tsv # 代谢通路覆盖度
# 将基因家族转换为EC号丰度
humann_regroup_table \
--input humann_output/sample_genefamilies.tsv \ # 输入：基因家族丰度
--output sample_ec.tsv \ # 输出：EC号丰度
--groups uniref90_level4ec # 按EC号重新分组
# 归一化
humann_renorm_table \
--input sample_ec.tsv \ # 输入
--output sample_ec_cpm.tsv \ # 输出
--units cpm # 归一化为CPM（每百万copies）
五、常见报错与解决¶
| 报错信息 | 原因 | 解决方案 |
|---|---|---|
dbCAN database not found | 数据库路径错误或未下载 | 重新运行dbcan_build或检查路径 |
HMMER: no hits | 蛋白序列太短或质量差 | 检查基因预测结果，过滤<100aa序列 |
DeepEC: CUDA out of memory | GPU显存不足 | 改用--device cpu或减小batch_size |
HUMAnN: database not found | 未下载HUMAnN数据库 | 运行humann_databases --download |
Prodigal: input error | 输入文件格式错误 | 确认是FASTA格式，检查是否有特殊字符 |
run_dbcan: permission denied | 数据库文件权限问题 | chmod -R 755 db/ |
六、面试高频问题¶
Q1: 从序列到酶活性预测的主要方法有哪些？¶
A: 三种主要方法：①基于同源性的方法（BLAST/DIAMOND比对已知酶数据库）；②基于HMM的方法（HMMER搜索酶家族特征模型）；③基于深度学习的方法（DeepEC/CLEAN直接从序列预测）。
Q2: CAZyme的六大类分别做什么？¶
A: GH(糖苷水解酶)=切糖链；GT(糖基转移酶)=接糖链；PL(多糖裂解酶)=裂解多糖；CE(碳水化合物酯酶)=去乙酰基等修饰；AA(辅助活性酶)=辅助降解(如木质素)；CBM(碳水化合物结合模块)=不是酶，帮助酶结合底物。
Q3: 宏基因组中预测的酶如何实验验证？¶
A: ①异源表达：克隆到大肠杆菌中表达并测活性；②宏转录组验证：确认基因在表达；③代谢组验证：检测预测的产物是否存在；④体外酶活测定。
七、速查表¶
# ===== CAZyme分析速查 =====
# 基因预测
prodigal -i assembly.fa -a proteins.faa -p meta
# dbCAN注释
run_dbcan proteins.faa protein --db_dir db/ --out_dir output/ --tools all
# 过滤标准：至少2种方法一致
# ===== EC号预测速查 =====
# DeepEC
python deepec.py -i proteins.faa -o output/
# HUMAnN3
humann --input reads.fq.gz --output output/ --threads 8
humann_regroup_table --input genefamilies.tsv --output ec.tsv --groups uniref90_level4ec
# ===== CAZyme大类速记 =====
# GH = 切糖 (Glycoside Hydrolase)
# GT = 接糖 (Glycosyl Transferase)
# PL = 裂糖 (Polysaccharide Lyase)
# CE = 去修饰 (Carbohydrate Esterase)
# AA = 辅助 (Auxiliary Activity)
# CBM = 粘住 (Carbohydrate-Binding Module)
# ===== EC号一级分类速记 =====
# 1 = 氧化还原 2 = 转移 3 = 水解
# 4 = 裂合 5 = 异构 6 = 连接 7 = 转位