KO丰度转换为pathway丰度
2020-05-01 本文已影响0人
kkkkkkang
经常会遇到已经有KO-number(代表基因家族信息)在各样品中的丰度，如何转换为level3等级的KEGG pathway(通路)丰度表呢？
目标如下：
KO丰度表
转换为
pathway丰度表
首先你要通过picrust或者metagenome一通计算得到第一个表，想得到第二个表有两种方法：
1. 借助picrust2，这个软件是通过16S 预测宏基因组功能的。很简单，就一行命令。注意--no_regroup参数必须加上，map的文件KEGG_pathways_to_KO.tsv找不到可以在家目录下find -name KEGG_pathways_to_KO.tsv
参考文献：https://www.frontiersin.org/articles/10.3389/fmicb.2019.01682/full
(picrust2) [jkyin@mn02 tables]$ pathway_pipeline.py -i PZH.KO.abund.copy.tsv -o kegg_pathway_out_abund --no_regroup --map ~/miniconda3/envs/picrust2/lib/python3.6/site-packages/picrust2/default_files/pathway_mapfiles/KEGG_pathways_to_KO.tsv
2. 借助humann2，这个软件可以在种水平对宏基因组和宏转录进行分析
- 首先你可以写脚本（样品比较多的时候，awk 把分隔符设置成tab然后按列输出就行了）或者简单的在excel中拆分成单个样本就行了，文件名一定要以.tsv结尾。把拆分后的样品文件放到一个文件夹中，当作humann2的输入（它可以接受多种格式文件的输入，比如raw reads，比对后的bam或者genetable等等），我们这里就属于genetable
- 比对需要humann1中的数据库，就直接把humann-0.99安装包下载下来解压就行了
(humann2) [jkyin@mn02 PZH_KEGG_humann2]$ for i in `ls`;do humann2 -i $i -o keggc/${i%%.*} --input-format genetable --pathways-database ~/biosoft/humann-0.99/data/keggc;done
- 然后把输出文件keggc/your/sample/*_pathabundance.tsv放到一个文件夹
(humann2) [jkyin@mn02 keggc]$ for i in `ls`; do cd $i;cp ${i}_pathabundance.tsv ../../pathabund/;cd ..;done
- 合并样品
(humann2) [jkyin@mn02 PZH_KEGG_humann2]$ humann2_join_tables -i pathabund/ -o all_join_pabd.tsv
Gene table created: /public/home/jkyin/yjk/PZH_meta/squeezemeta/PZH/results/tables/PZH_KEGG_humann2/all_join_pabd.tsv
好，现在两种方法殊途同归。humann2多了UNMAPPED和UNINTEGRATED，没关系，反正我们也不care这俩。humann2得到的pathway比picrust2多几个，对结果基本没有影响。
picrust2结果
humann2结果
只有ko00010这种pathway编号，没有名字怎么办呢？humann2来帮你
- 重命名以及cpm（counts per million）标准化，或者你也可以进行relative abundance 标准化
(humann2) [jkyin@mn02 PZH_KEGG_humann2]$ humann2_rename_table -i all_join_pabd.tsv -o all_join_pabd_rename.tsv -n kegg-pathway
Loading table from: all_join_pabd.tsv
Loading mapping file from: /public/home/jkyin/miniconda3/envs/humann2/lib/python2.7/site-packages/humann2/tools/../data/misc/map_kegg-pwy_name.txt.gz
Renamed 251 of 253 entries (99.21%)
(humann2) [jkyin@mn02 PZH_KEGG_humann2]$ humann2_renorm_table -i all_join_pabd_rename.tsv -o all_join_pabd_rename_cpm.tsv -u cpm Loading table from: all_join_pabd_rename.tsv
好啦，完成~
