1. 这个工具有在线网页版，当然KEGG注释还有其他很多工具可以选择比如Diamond，KAAS等等
https://www.genome.jp/tools/kofamkoala/
2. 安装和使用都可以参考，这里主要是后续分析
https://blog.csdn.net/woodcorpse/article/details/106554548
最后会输出
输出文件，--cpu给了20，内存给了40G，跑了块4天
到这思路就很简单了，因为我们已经拿到了各个基因在每个样本中的相对丰度，转化一下就是这些K号在不同样本中的相对丰度，然后加起来即可。然后K号转为level1 level2 level3
cat Prodigal_KEGG_raw.txt | tr "\t" "," > Prodigal_KEGG_raw.csv 先转成逗号分隔符
如果用genemark的结果会报错，但是检查过，这个名字的的确确是unique的，很怪，待解决
3. 获取KO号的相对丰度表
随便写了个脚本
python KOAbdun.py Prodigal_KEGG_raw.csv > KO_Abun.raw.csv #Prodigal_KEGG_raw.csv即KofamKOALA的输出文件，不过我这里转成了逗号分隔符
~~~
awk 'BEGIN{ FS=",";OFS="," }{a[$1]+=$2}{b[$1]+=$3}{c[$1]+=$4}{d[$1]+=$5}{e[$1]+=$6}{f[$1]+=$7}{g[$1]+=$8}{h[$1]+=$9}{m[$1]+=$10}{j[$1]+=$11}{k[$1]+=$12}{l[$1]+=$13}END{for(i in a){print i,a[i],b[i],c[i],d[i],e[i],f[i],g[i],h[i],m[i],j[i],k[i],l[i]}}' KO_Abun.raw.csv | sort -k 1 > tmp
# awk的作用是合并，如果第一列相同，就把后面12项相加，因为我有12个样本，合并为一行
head -n 1 KO_Abun.raw.csv > head ; cat head tmp > KO_Abun.csv
~~~
KO_Abun.csv即最终的KO号的相对丰度表
4. KEGG Level1/2/3的注释
再看下Kofam输出文件
把空行删掉
4.1我在网上找到了KO号和pathway_id的对应关系 (这个真是找麻了，有需要可以找我要
4.2又找到了pathway_id和Level1/2/3的对应关系 Ref: https://www.jianshu.com/p/beb0c3727a15
把文件处理成这样：
kegg_pathway.txt
因为每个基因的相对丰度都是有的，所以map号的相对丰度也就有了
比如对于功能A，有两个基因注释到了，那功能A的相对丰度就是这两个基因的相对丰度之和
python gene2K2mapID.py kegg_pathway.txt > map_Abun.csv
~~~
因为有不少gene注释到了同一个map号
所以用awk合并相同行
awk 'BEGIN{ FS="\t";OFS="\t" }{a[$1]+=$2}{b[$1]+=$3}{c[$1]+=$4}{d[$1]+=$5}{e[$1]+=$6}{f[$1]+=$7}{g[$1]+=$8}{h[$1]+=$9}{m[$1]+=$10}{j[$1]+=$11}{k[$1]+=$12}{l[$1]+=$13}END{for(i in a){print i,a[i],b[i],c[i],d[i],e[i],f[i],g[i],h[i],m[i],j[i],k[i],l[i]}}' map_Abun.csv > map_Abun.csv3
这行代码的意思是 如果第一列相同，则其余列相加，因为我有12个样，所以除了第一列，后面有12列数据，所以这里字母有12个 a[i] 到 l[i]
~~~
python mapID2Lv123.py map_Abun.csv3 > Lv123.Abun.tsv2
第一列是map号，第二列是Level1，第三列是Level2，第四列是Level4，到这其实就很简单了
5.最后拿到了Level1/2/3的相对丰度表，随便画个堆叠柱状图
用LEfSe找差异功能： https://www.jianshu.com/p/35e3f725c554
5.5关于KO号的可视化
只能找到这些了，一个α多样性，一个β多样性的降维，一个热图，就是这个热图咋画的我还不是很明白
Xiong et al. Microbiome (2021) 9:171
https://doi.org/10.1186/s40168-021-01118-6
《Plant developmental stage drives the differentiation in ecological role of the maize microbiome》
2022.3更新
由于K号(eg:K00174)和pathway id之间并不是一一对应的关系，而是被包含和包含 这样一个子集父集的关系，所以简单粗暴的把K号和pathway id进行一 一对应转换，不够合理
发现了一个ReportScore的工具，可以解决问题，正是利用被包含和包含的关系，类比于功能富集，不过这个时候的背景就是KEGG pathway本身，输入的基因也就是我们得到的K号
详见：
https://mp.weixin.qq.com/s/3RCxhQk357D7PShfRXTINg
https://github.com/wangpeng407/ReporterScore
6.脚本
gene2K2mapID.py
mapID2Lv123.py
pathway_annotation.txt即map号和Level123的对应关系
.
网友评论