導(dǎo)讀

1 MetaPhlAn獲得宏基因組物種和物種豐度
2 MetaPhlAn結(jié)果合并和可視化
3 GraPhlAn繪制進(jìn)化樹
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數(shù)據(jù)

cd metaphlan
mkdir input

# 10個(gè)口腔粘膜樣品
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個(gè)舌背樣品
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. 結(jié)果

profiled_samples文件夾生成20個(gè)結(jié)果文件，每個(gè)文件格式如下：

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結(jié)果合并和可視化

依賴：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 標(biāo)準(zhǔn)化方法
--top 25 default --perc 90 物種選取方案
-c bbcry 配色方案
default -m average 聚類方案
default -d braycurtis 舉例計(jì)算方法
default -f correlation 計(jì)算樣品相關(guān)

三、GraPhlAn可視化

1. 依賴：graphlan matplotlib

hg clone ssh://hg@bitbucket.org/nsegata/graphlan

2. 單樣品準(zhǔn)備和可視化

# 準(zhǔn)備
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. 多樣品準(zhǔn)備和可視化

# 準(zhǔn)備
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的結(jié)果導(dǎo)入krona進(jìn)行可視化
krona: http://sourceforge.net/p/krona/home/krona/

四、LEfSe尋找分類biomarker

4.1 數(shù)據(jù)準(zhǔn)備

# 1. 數(shù)據(jù)整理
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 繪制進(jìn)化樹圖

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打包結(jié)果】

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蛋白庫(kù)

# 獲取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