PyKofamSearch：基于 PyHMMER 的快速实现KEGG KO 功能注释

# 创建一个名为PyKofamSearch的conda环境，并安装pykofamsearch
conda create -n PyKofamSearch -c conda-forge -c bioconda pykofamsearch
# conda install 直接安装
conda install -c bioconda pykofamsearch
# pip install 安装
pip install pykofamsearch
# 激活环境
conda activate PyKofamSearch  

## 下载数据库
# 在线模式
serialize_kofam_models -o db/KEGG
# 离线模式
mkdir -p db/KEGG && cd db/KEGG 
# 使用wget下kofam数据库Hmm文件和ko注释文件
wget -c -nv -a download.log -b https://www.genome.jp/ftp/db/kofam/profiles.tar.gz 
wget -c -nv -a download.log -b https://www.genome.jp/ftp/db/kofam/ko_list.gz 
# 解压数据库
gzip -d ko_list.gz 
tar -xvzf profiles.tar.gz 
# 格式化数据库
serialize_kofam_models -d profiles -k ko_list.tsv -b KofamSearch.pkl.gz

pykofamsearch -h
usage: pykofamsearch -i <proteins.fasta> -o <output.tsv> -d 

    Running: pykofamsearch v2025.9.5 via Python v3.12.13 | /path/softwara/miniforge3/envs/PyKofamSearch/bin/python3.12

options:
  -h, --help            show this help message and exit
  --verbosity VERBOSITY
                        Verbosity of missing KOfams [Default: 1]
  -v, --version         show program's version number and exit

I/O arguments:
  -i PROTEINS, --proteins PROTEINS
                        path/to/proteins.fasta. stdin does not stream and loads everything into memory. [Default: stdin]
  -o OUTPUT, --output OUTPUT
                        path/to/output.tsv [Default: stdout]
  -s SUBSET, --subset SUBSET
                        path/to/identifiers.list where HMM identifiers are on a separate line used to subset the database. Only HMMs in the subset will be used.
  --no_header           No header

Utility arguments:
  -p N_JOBS, --n_jobs N_JOBS
                        Number of threads to use [Default: 1]

HMMSearch arguments:
  -e EVALUE, --evalue EVALUE
                        E-value threshold [Default: 0.1]
  -a, --all_hits        Return all hits and do not use curated threshold. Not recommended for large queries.
  -t THRESHOLD_SCALE, --threshold_scale THRESHOLD_SCALE
                        Multiplier for the curated thresholds. Higher values will make the annotation more strict [Default: 1.0]

Database arguments:
  -d DATABASE_DIRECTORY, --database_directory DATABASE_DIRECTORY
                        path/to/kofam_database_directory/ cannot be used with -b/-serialized_database
  -b SERIALIZED_DATABASE, --serialized_database SERIALIZED_DATABASE
                        path/to/database.pkl cannot be used with -d/--database_directory

PyKOfamSearch

# 示例蛋白序列：https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_protein.faa.gz
# 去除序列ID后的注释信息
seqkit replace -p '\s.*' -r '' GCF_000005845.2_ASM584v2_protein.faa.gz -w 0 > ASM584v2_protein.faa
# 使用默认参数对蛋白进行注释, 默认依据 KOfam 数据库为各 KO 模型设定的 curated threshold 对命中结果进行筛选，仅保留 score 达到或超过对应 threshold 的显著注释结果。
pykofamsearch -i ASM584v2_protein.faa -b ./db/KEGG/PyKofamSearch.pkl.gz -o ASM584v2_kegg.tsv -p 16
# 使用 -a 参数返回所有满足 E-value 阈值的 HMM 命中结果将返回所有满足 E-value 阈值的 HMM 命中结果，不再依据 KOfam curated threshold 对结果进行筛选，因此输出中会同时包含 score 高于和低于对应 threshold 的命中结果。
pykofamsearch -i ASM584v2_protein.faa -b ./db/KEGG/PyKofamSearch.pkl.gz -o ASM584v2_kegg.tsv -p 16 -a
# 使用 -a 参数返回所有满足 E-value 阈值的 HMM 命中结果，并通过 -e 0.001 将 E-value 阈值设为 0.001，即保留 E-value ≤ 0.001 的命中结果，不再依据 KOfam curated threshold 进行筛选。
pykofamsearch -i ASM584v2_protein.faa -b ./db/KEGG/PyKofamSearch.pkl.gz -o ASM584v2_kegg.tsv -p 16 -a -e 0.001

# 查看结果
csvtk pretty -t ASM584v2_kegg.tsv | head