宏基因组 - (3)使用KofamKOALA对非冗余基因集进行K

2021-12-17 本文已影响0人深山夕照深秋雨OvO

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