记录细菌基因组KEGG注释方法。kofanscan可源码安装（参考链接）或者conda安装（无需管理员权限，推荐）。数据库需要额外单独下载（参考链接）。
文章
标题：KofamKOALA: KEGG ortholog assignment based on profile HMM and adaptive score threshold
期刊：Bioinformatics
时间：2020
被引：239（谷歌学术 2021.11.20）
prokka获取细菌基因组蛋白预测文件
conda activate /route/miniconda3/envs/prokka
prokka 22607.fasta \
--outdir prokka \
--prefix 22607 \
--cpus 16 \
--kingdom Bacteria --force
kofamkoala注释KEGG
官网：KofamKOALA - KEGG Orthology Search
数据库：https://www.genome.jp/ftp/db/kofam/
下载、安装：KEGG功能注释工具 KofamKOALA 安装与使用
安装依赖，再安装主程序
# 安装依赖，再安装主程序
conda install ruby
conda install parallel
conda install hmmer
conda install kofamscan
exec_annotation --help
直接安装也可以
conda create -n kofam
conda activate kofam
conda install kofamscan -c bioconda
# 安装过程：
kofamscan version 1.3.0-2 has been successfully installed!
This software needs a database which can be downloaded from ftp.genome.jp/pub/db/kofam/
For more details see ftp://ftp.genome.jp/pub/tools/kofamscan/README
done
KEGG注释
source /hwfssz5/ST_META/P18Z10200N0423_ZYQ/MiceGutProject/hutongyuan/hwfssz1/software/miniconda3/etc/profile.d/conda.sh
conda activate r403
mkdir ko_tmp
exec_annotation \
#-f mapper \
-f detail-tsv \
-E 1e-3 \
--profile /route/hutongyuan/database/kofamscan/profiles \
--ko-list /route/hutongyuan/database/kofamscan/ko_list \
--cpu 6 \
--tmp-dir ko_tmp \
-o ./kegg_detail_neg3.txt \
../prokka/22607.faa
rm -r ko_tmp
结果查看
NMONGKEM_00013 K07574
NMONGKEM_00014 K06948
NMONGKEM_00015 K07015
NMONGKEM_00016 K02887
NMONGKEM_00017 K02916
NMONGKEM_00018 K02520
NMONGKEM_00019 K01868
...
-f mapper：结果输出只保留query蛋白序列id和ko id
-f detail-tsv：完全比对结果
EC信息提取
# 提取全部EC
for i in `ls kegg_detail/`; do
cat kegg_detail/$i | grep 'EC:' | awk -F"EC:" '{print $2}' | sed 's/]"//' | grep -v '-' | sed 's/ /\n/g' | sort | uniq > kegg_ec_detail/$i
echo -e "$i done..."
done
# 提取score > threshold的EC
for i in `ls kegg_detail/`; do
cat kegg_detail/$i | grep 'EC:' | awk -F"\t" '{if($5>$4) print $0}' | awk -F"EC:" '{print $2}' | sed 's/]"//' | grep -v '-' | sed 's/ /\n/g' | sort | uniq > kegg_ec_high/$i
echo -e "$i done..."
done
# grep 查找有EC的行
# awk 判断score是否大于阈值
# awk 拆分取出EC信息
# sed 把格式符去掉
# grep -v 去掉-行
# sed 一对多换行
# sort uniq 排序，去重
二、注释module
利用ko2module对应关系进行匹配
K00844 M00001
K12407 M00001
K00845 M00001
K00886 M00001
K08074 M00001
三、注释pathway
利用ko2pathway对应关系进行匹配
K00001 map00010
K00002 map00010
K00016 map00010
K00114 map00010
K00121 map00010
K00128 map00010
四、分层注释pathway
利用map_id2level对应关系匹配
01100 Metabolism Global and overview maps Metabolic pathways
01110 Metabolism Global and overview maps Biosynthesis of secondary metabolites
01120 Metabolism Global and overview maps Microbial metabolism in diverse environments
01200 Metabolism Global and overview maps Carbon metabolism
比较humann3和kofam的数据库大小
# humann3 2021.9
zcat map_ko_name.txt.gz | wc -l
23326
# kofam 2021.9
cat ko_list | wc -l
24541
差不多大，kofam险胜
Module功能模块，是对KO的分类，涉及多种通路map和反应
更多
一个非模式生物的KEGG注释和使用clusterProfiler做KEGG富集的实例
library(clusterProfiler)
setwd('G:/genome-Cal/write/CalvsEla/KEGG/')
kegg = read.csv('kegg.anno.xls',header = T,sep = '\t')
gene_list = read.csv('gene.list',header = F,sep = '\t')
gene_list = gene_list$V1
enrich_KEGG = enricher(gene = gene_list,TERM2GENE = kegg[c('ko_id','gene_id')],TERM2NAME = kegg[c('ko_id','pathway')],pvalueCutoff = 0.5)
dotplot(enrich_KEGG)
#kegg.anno.xls 上面得到的KEGG背景文件
#gene.list 要富集的gene id列表，如差异基因
