导读
1 MetaPhlAn获得宏基因组物种和物种丰度
2 MetaPhlAn结果合并和可视化
3 GraPhlAn绘制进化树
4 LEfSe寻找分类biomarker并可视化
5 HUMAnN代谢分析
bitbucket地址：https://bitbucket.org/nsegata/metaphlan/wiki/MetaPhlAn_Pipelines_Tutorial
一、MetaPhlAn分析 20 HMP 样品
1. 下载安装MetaPhlAn
wget https://bitbucket.org/nsegata/metaphlan/get/default.tar.bz2
tar xjvf default.tar.bz2
mv *-metaphlan-* metaphlan
2. 获取20 HMP数据
cd metaphlan
mkdir input
# 10个口腔粘膜样品
buccal_mucosa_samples="SRS013506 SRS015374 SRS015646 SRS017687 SRS019221 SRS019329 SRS020336 SRS022145 SRS022532 SRS045049"
for s in ${buccal_mucosa_samples}
do
wget http://downloads.hmpdacc.org/data/Illumina/buccal_mucosa/${s}.tar.bz2 -O input/${s}.tar.bz2
done
# 10个舌背样品
tongue_dorsum_samples="SRS011243 SRS013234 SRS014888 SRS015941 SRS016086 SRS016342 SRS017713 SRS019219 SRS019327 SRS043663"
for s in ${tongue_dorsum_samples}
do
wget http://downloads.hmpdacc.org/data/Illumina/tongue_dorsum/${s}.tar.bz2 -O input/${s}.tar.bz2
done
3. MetaPhlAn分析
mkdir profiled_samples
# 口腔粘膜样品
BM_samples="SRS013506 SRS015374 SRS015646 SRS017687 SRS019221 SRS019329 SRS020336 SRS022145 SRS022532 SRS045049"
for s in ${BM_samples}
do
tar xjf input/${s}.tar.bz2 --to-stdout | ./metaphlan.py --bowtie2db bowtie2db/mpa --bt2_ps very-sensitive --input_type multifastq --bowtie2out BM_${s}.bt2out > profiled_samples/BM_${s}.txt
done
# 舌背样品
TD_samples="SRS011243 SRS013234 SRS014888 SRS015941 SRS016086 SRS016342 SRS017713 SRS019219 SRS019327 SRS043663"
for s in ${TD_samples}
do
tar xjf input/${s}.tar.bz2 --to-stdout | ./metaphlan.py --bowtie2db bowtie2db/mpa --bt2_ps very-sensitive --input_type multifastq --bowtie2out TD_${s}.bt2out > profiled_samples/TD_${s}.txt
done
4. 结果
profiled_samples文件夹生成20个结果文件，每个文件格式如下：
cat profiled_samples/BM_SRS013506.txt
$ k__Bacteria 100.0
$ k__Bacteria|p__Firmicutes 80.97874
$ k__Bacteria|p__Proteobacteria 17.17081
$ k__Bacteria|p__Fusobacteria 0.34233
$ k__Bacteria|p__Bacteroidetes 0.17203
$ k__Bacteria|p__Firmicutes|c__Bacilli 80.62406
[truncated output]
二、MetaPhlAn结果合并和可视化
依赖：matplotlib python包
2.1 合并文件
mkdir output
utils/merge_metaphlan_tables.py profiled_samples/*.txt > output/merged_abundance_table.txt
2.2 可视化：heatmap图
mkdir output_images
plotting_scripts/metaphlan_hclust_heatmap.py \
-c bbcry --top 25 --minv 0.1 -s log \
--in output/merged_abundance_table.txt \
--out output_images/abundance_heatmap.png
default --tax_lev s 仅展示物种
-s log 标准化方法
--top 25 default --perc 90 物种选取方案
-c bbcry 配色方案
default -m average 聚类方案
default -d braycurtis 举例计算方法
default -f correlation 计算样品相关
三、GraPhlAn可视化
1. 依赖：graphlan matplotlib
hg clone ssh://hg@bitbucket.org/nsegata/graphlan
2. 单样品准备和可视化
# 准备
mkdir tmp
plotting_scripts/metaphlan2graphlan.py profiled_samples/BM_SRS013506.txt --tree_file tmp/BM_SRS013506.tree.txt --annot_file tmp/BM_SRS013506.annot.txt
# 可视化
graphlan_annotate.py --annot tmp/BM_SRS013506.annot.txt tmp/BM_SRS013506.tree.txt tmp/BM_SRS013506.xml
graphlan.py --dpi 200 tmp/BM_SRS013506.xml output_images/BM_SRS013506.png
2. 多样品准备和可视化
# 准备
mkdir -p tmp
for file in profiled_samples/*
do
filename=`basename ${file}`
samplename=${filename%\.*}
plotting_scripts/metaphlan2graphlan.py ${file} --tree_file tmp/${samplename}.tree.txt --annot_file tmp/${samplename}.annot.txt
graphlan_annotate.py --annot tmp/${samplename}.annot.txt tmp/${samplename}.tree.txt tmp/${samplename}.xml
graphlan.py --dpi 200 tmp/${samplename}.xml output_images/${samplename}.png
done
# 可视化
plotting_scripts/metaphlan2graphlan.py output/merged_abundance_table.txt --tree_file tmp/merged.tree.txt --annot_file tmp/merged.annot.txt
graphlan_annotate.py --annot tmp/merged.annot.txt tmp/merged.tree.txt tmp/merged.xml
graphlan.py --dpi 200 tmp/merged.xml output_images/merged.png
另外，使用conversion_scripts文件夹中的metaphlan2krona.py脚本可以直接将metaphlan的结果导入krona进行可视化
krona: http://sourceforge.net/p/krona/home/krona/
四、LEfSe寻找分类biomarker
4.1 数据准备
# 1. 数据整理
sed 's/\([A-Z][A-Z]\)_\w*/\1/g' output/merged_abundance_table.txt > tmp/merged_abundance_table.4lefse.txt
# 2. lefse格式化
format_input.py tmp/merged_abundance_table.4lefse.txt tmp/merged_abundance_table.lefse -c 1 -o 1000000
# 3. lefse分析
run_lefse.py tmp/merged_abundance_table.lefse tmp/merged_abundance_table.lefse.out -l 4
4.2 绘制LDA图
plot_res.py --dpi 300 tmp/merged_abundance_table.lefse.out output_images/lefse_biomarkers.png
4.3 绘制进化树图
plot_cladogram.py --dpi 300 --format png tmp/merged_abundance_table.lefse.out output_images/lefse_biomarkers_cladogram.png
4.4 单物种组间差异柱形图
plot_features.py -f one --feature_name "k__Bacteria.p__Firmicutes" tmp/merged_abundance_table.lefse tmp/merged_abundance_table.lefse.out output/Firmicutes.png
4.5 所有物种组间差异柱形图【zip打包结果】
plot_features.py -f diff --archive zip tmp/merged_abundance_table.lefse tmp/merged_abundance_table.lefse.out biomarkers.zip
五、HUMAnN代谢分析
1. 依赖：HUMAnN、scon、KEGG蛋白库
# 获取HUMAnN
hg clone ssh://hg@bitbucket.org/chuttenh/humann
2. HUMAnN使用
mkdir output
utils/merge_metaphlan_tables.py humann_profiling/*.txt > output/merged_humann_abundance_table.txt
