MetaPhlAn 4.0 tutorial
MetaPhlAn4是一个从宏基因组shotgun测序数据中分析微生物群落组成的工具
数据库位置 Folder containing the MetaPhlAn database. You can specify the location by exporting the DEFAULT_DB_FOLDER variable in the shell.
--bowtie2db METAPHLAN_BOWTIE2_DB [default /xxx/.conda/envs/metagenome/lib/python3.7/site-packages/metaphlan/metaphlan_databases]
MetaPhlAn的输入文件可以是fastq文件：SRS014476-Supragingival_plaque.fasta.gz
>HWUSI-EAS1568_102539179:1:100:10001:9884/1 TTCCTTGCTGTAATAGGTGCCACGATGGTCGGATTGATATTGCCAACTGTTATCTGTTCGCACAATGATAAAGCCTTGTTTCGGTTGTGGCGG >HWUSI-EAS1568_102539179:1:100:10002:15186/1 CCTTGCTCATCTTTTACCAAATAAACCACGGCCTGTTTCGCTCGCACTTTTAAGCCTGCTTTATCCTGGTCAGCAGGAATCGGCTGACTAAG
运行metaphlan
metaphlan SRS014476-Supragingival_plaque.fasta.gz --input_type fasta > SRS014476-Supragingival_plaque_profile.txt
运行结果SRS014476-Supragingival_plaque.fasta.gz.bowtie2out.txt
SRS014476-Supragingival_plaque.fasta.gz.bowtie2out.txt
SRS014476-Supragingival_plaque_profile.txt
#mpa_vJan21_CHOCOPhlAnSGB_202103表示 MetaPhlAn 在这次运行中使用的标记基因数据库的版本。目前在数据库中包括从～100k 参考基因组(～99,500个细菌和古细菌和? 500个真核生物)中鉴定的～ 1.1 M 独特的进化枝特异性标记基因。clade_name列：这一列报道的分类谱系包括界(例如细菌/古细菌)和物种基因组仓(SGB)。分类名称的前缀有助于表明它们的等级Kingdom: k__, Phylum: p__, Class: c__, Order: o__, Family: f__, Genus: g__, Species: s__, and SGB: t__.relative_abundance(%)。由于典型的shutgun测序为基础的分类概况是相对的(即它不提供绝对细胞计数) ，进化枝是分层总结。每个分类层次的总和将达到100% 。也就是说，所有界级别的进化枝的总和是100% ，所有门级别的进化枝的总和是100% ，以此类推。
#mpa_vJan21_CHOCOPhlAnSGB_202103
clade_name
relative_abundance
查看所有目级别的分类总和是100%
grep o__ SRS014476-Supragingival_plaque_profile.txt | grep -v f__ | cut -f1,3
gitee地址toy_databases
databases/ ├── mpa_latest ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.1.bt2l ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.2.bt2l ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.3.bt2l ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.4.bt2l ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.fna ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.pkl ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.rev.1.bt2l ├── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.rev.2.bt2l └── mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.tar
metaphlan \ --bowtie2db databases \ --index mpa_vJan21_TOY_CHOCOPhlAnSGB_202103 \ SRS014476-Supragingival_plaque.fasta.gz --input_type fasta > SRS014476-Supragingival_plaque_profile.txt
WARNING: The number of reads in the sample (None) is below the recommended minimum of 10,000 reads.WARNING: MetaPhlAn did not detect any microbial taxa in the sample.
小技巧：循环批量处理样本列表
# 基于样本元数据提取样本列表命令解析 # 去掉表头 tail -n+2 result/metadata.txt # 提取第一列样本名 tail -n+2 result/metadata.txt|cut -f1 # 循环处理样本 for i in `tail -n+2 result/metadata.txt|cut -f1`;do echo "Processing "$i; done # ` 反引号为键盘左上角Esc键下面的按键，一般在数字1的左边，代表运行命令返回结果
HUMAnN2要求双端序列合并的文件作为输入，for循环根据实验设计样本名批量双端序列合并。注意星号和问号，分别代表多个和单个字符。当然大家更不能溜号~~~行分割的代码行末有一个 \
mkdir -p temp/concat # 双端合并为单个文件 for i in `tail -n+2 result/metadata.txt|cut -f1`;do cat temp/qc/${i}_?.fastq \ > temp/concat/${i}.fq; done # 查看样品数量和大小 ls -shl temp/concat/*.fq # 数据太大，计算时间长，可用head对单端分析截取20M序列，即3G，则为80M行，详见附录：HUMAnN2减少输出文件加速
启动humann2环境：仅humann2布置于自定义环境下使用
# 方法1. conda加载环境 conda activate humann2 # 方法2. source加载指定 # source /home/liuyongxin/miniconda2/envs/humann2/bin/activate
检查数据库配置是否正确
humann2 --version # v2.8.1 humann2_config mkdir -p temp/humann2
单样本1.25M PE150运行测试，8p，2.5M，1~2h；0.2M, 34m；0.1M，30m；0.01M，25m；16p，18m
# CRITICAL ERROR: Can not call software version for bowtie2，见"Perl环境" i=C1 # memusg -t humann2 --input temp/concat/${i}.fq --output temp/humann2 --threads 16
多样本并行计算，测试数据约30m，系统耗时12小时
tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \ 'humann2 --input temp/concat/{1}.fq \ --output temp/humann2/' # (可选)大文件清理，humann2临时文件可达原始数据30~40倍 # 链接重要文件至humann2目录 for i in `tail -n+2 result/metadata.txt|cut -f1`;do ln temp/humann2/${i}_humann2_temp/${i}_metaphlan_bugs_list.tsv temp/humann2/ done # 删除临时文件 rm -rf temp/concat/* temp/humann2/*_humann2_temp
mkdir -p result/metaphlan2 # 合并、修正样本名、预览 merge_metaphlan_tables.py temp/humann2/*_metaphlan_bugs_list.tsv | \ sed 's/_metaphlan_bugs_list//g' > result/metaphlan2/taxonomy.tsv head -n5 result/metaphlan2/taxonomy.tsv
metaphlan_to_stamp.pl result/metaphlan2/taxonomy.tsv \ > result/metaphlan2/taxonomy.spf head -n5 result/metaphlan2/taxonomy.spf # 下载metadata.txt和taxonomy.spf使用stamp分析 # 网络分析见附录 metaphlan2-共有或特有物种网络图
# c设置颜色方案，top设置物种数量，minv最小相对丰度，s标准化方法，log为取10为底对数，xy为势图宽和高，图片可选pdf/png/svg格式 metaphlan_hclust_heatmap.py \ --in result/metaphlan2/taxonomy.tsv \ --out result/metaphlan2/heatmap.pdf \ -c jet --top 30 --minv 0.1 \ -s log -x 0.4 -y 0.2 # 报错解决详见附录：### metaphlan_hclust_heatmap.py报错AttributeError: Unknown property axisbg # 帮助见 metaphlan_hclust_heatmap.py -h # 更多绘制见3StatPlot.sh
合并通路丰度(pathabundance)，含功能和对应物种组成。可选基因家族(genefamilies 太多)，通路覆盖度(pathcoverage)。注意看屏幕输出# Gene table created: result/humann2/pathabundance.tsv
# Gene table created: result/humann2/pathabundance.tsv
mkdir -p result/humann2 humann2_join_tables \ --input temp/humann2 \ --file_name pathabundance \ --output result/humann2/pathabundance.tsv # 样本名调整：删除列名多余信息 head result/humann2/pathabundance.tsv sed -i 's/_Abundance//g' result/humann2/pathabundance.tsv # 预览和统计 head result/humann2/pathabundance.tsv csvtk -t stat result/humann2/pathabundance.tsv
标准化为相对丰度relab(1)或百万比cpm(1,000,000)
humann2_renorm_table \ --input result/humann2/pathabundance.tsv \ --units relab \ --output result/humann2/pathabundance_relab.tsv head -n5 result/humann2/pathabundance_relab.tsv
分层结果：功能和对应物种表(stratified)和功能组成表(unstratified)
humann2_split_stratified_table \ --input result/humann2/pathabundance_relab.tsv \ --output result/humann2/ # 可以使用stamp进行统计分析
两样本无法组间比较，在pcl层面替换为HMP数据进行统计和可视化。
参考 https://bitbucket.org/biobakery/humann2/wiki/Home#markdown-header-standard-workflow
在通路丰度中添加分组
## 提取样品列表 head -n1 result/humann2/pathabundance.tsv | sed 's/# Pathway/SampleID/' | tr '\t' '\n' > temp/header ## 对应分组，本示例分组为第2列($2)，根据实际情况修改 awk 'BEGIN{FS=OFS="\t"}NR==FNR{a[$1]=$2}NR>FNR{print a[$1]}' result/metadata.txt temp/header | tr '\n' '\t'|sed 's/\t$/\n/' > temp/group # 合成样本、分组+数据 cat <(head -n1 result/humann2/pathabundance.tsv) temp/group <(tail -n+2 result/humann2/pathabundance.tsv) \ > result/humann2/pathabundance.pcl head -n5 result/humann2/pathabundance.pcl
组间比较，样本量少无差异，结果为4列的文件：通路名字，通路在各个分组的丰度，差异P-value，校正后的Q-value。演示数据2样本无法统计，此处替换为HMP的结果演示统计和绘图(上传hmp_pathabund.pcl，替换pathabundance.pcl为hmp_pathabund.pcl)。
wget http://210.75.224.110/db/train/meta/result/humann2/hmp_pathabund.pcl # cp /db/humann2/hmp_pathabund.pcl ./ mv hmp_pathabund.pcl result/humann2/ # 设置输入文件名 pcl=result/humann2/hmp_pathabund.pcl # 统计表格行、列数量 csvtk -t stat ${pcl} head -n2 ${pcl} |cut -f 1-5 # 按分组KW检验，注意第二列的分组列名 humann2_associate --input ${pcl} \ --focal-metadatum Group --focal-type categorical \ --last-metadatum Group --fdr 0.05 \ --output result/humann2/associate.txt wc -l result/humann2/associate.txt head -n5 result/humann2/associate.txt
barplot展示通路的物种组成，如：腺苷核苷酸合成
# --sort sum metadata 按丰度和分组排序 # 指定差异通路，如 P163-PWY / PWY-3781 / PWY66-409 / PWY1F-823 path=PWY-3781 humann2_barplot --sort sum metadata \ --input ${pcl} --focal-feature ${path} \ --focal-metadatum Group --last-metadatum Group \ --output result/humann2/barplot_${path}.pdf
需要下载utility_mapping数据库并配置成功才可以使用。详见软件和数据库安装1soft_db.sh。
支持GO、PFAM、eggNOG、level4ec、KEGG的D级KO等注释，详见humann2_regroup_table -h。
humann2_regroup_table -h
# 转换基因家族为KO(uniref90_ko)，可选eggNOG(uniref90_eggnog)或酶(uniref90_level4ec) for i in `tail -n+2 result/metadata.txt|cut -f1`;do humann2_regroup_table \ -i temp/humann2/${i}_genefamilies.tsv \ -g uniref90_ko \ -o temp/humann2/${i}_ko.tsv done # 合并，并修正样本名 humann2_join_tables \ --input temp/humann2/ \ --file_name ko \ --output result/humann2/ko.tsv sed -i '1s/_Abundance-RPKs//g' result/humann2/ko.tsv tail result/humann2/ko.tsv # 与pathabundance类似，可进行标准化renorm、分层stratified、柱状图barplot等操作
# metaphlan2 to graphlan export2graphlan.py --skip_rows 1,2 -i result/metaphlan2/taxonomy.tsv \ --tree temp/merged_abundance.tree.txt \ --annotation temp/merged_abundance.annot.txt \ --most_abundant 1000 --abundance_threshold 20 --least_biomarkers 10 \ --annotations 3,4 --external_annotations 7 # 参数说明见PPT，或运行 export2graphlan.py --help # graphlan annotation graphlan_annotate.py --annot temp/merged_abundance.annot.txt \ temp/merged_abundance.tree.txt temp/merged_abundance.xml # output PDF figure, annoat and legend graphlan.py temp/merged_abundance.xml result/metaphlan2/graphlan.pdf \ --external_legends
前面演示数据仅有2个样本，无法进行差异比较。下面使用result12目录中由12个样本生成的结果表进行演示
# 设置结果目录，自己的数据使用result，演示用result12 result=result12 # 下载演示数据 # wget http://210.75.224.110/db/EasyMetagenome/result12.zip && unzip result12.zip
准备输入文件，修改样本品为组名(可手动修改)
# 预览输出数据 head -n3 $result/metaphlan2/taxonomy.tsv # 提取样本行，替换为每个样本一行，修改ID为SampleID head -n1 $result/metaphlan2/taxonomy.tsv|tr '\t' '\n'|sed '1 s/ID/SampleID/' >temp/sampleid head -n3 temp/sampleid # 提取SampleID对应的分组Group(假设为metadata.txt中第二列$2)，替换换行\n为制表符\t，再把行末制表符\t替换回换行 awk 'BEGIN{OFS=FS="\t"}NR==FNR{a[$1]=$2}NR>FNR{print a[$1]}' $result/metadata.txt temp/sampleid|tr '\n' '\t'|sed 's/\t$/\n/' >groupid cat groupid # 合并分组和数据(替换表头) cat groupid <(tail -n+2 $result/metaphlan2/taxonomy.tsv) > $result/metaphlan2/lefse.txt head -n3 $result/metaphlan2/lefse.txt
方法1. 推荐在线 http://www.ehbio.com/ImageGP 中LEfSe一键分析
方法2. (可选)LEfSe命令行分析代码供参考
# 格式转换为lefse内部格式 lefse-format_input.py \ $result/metaphlan2/lefse.txt \ temp/input.in -c 1 -o 1000000 # 运行lefse(样本无重复、分组将报错) run_lefse.py temp/input.in \ temp/input.res # 绘制物种树注释差异 lefse-plot_cladogram.py temp/input.res \ result/metaphlan2/lefse_cladogram.pdf --format pdf # 绘制所有差异features柱状图 lefse-plot_res.py temp/input.res \ $result/metaphlan2/lefse_res.pdf --format pdf # 绘制单个features柱状图 # 查看显著差异features，按丰度排序 grep -v '-' temp/input.res | sort -k3,3n # 手动选择指定feature绘图，如Firmicutes lefse-plot_features.py -f one \ --feature_name "k__Bacteria.p__Firmicutes" \ --format pdf \ temp/input.in temp/input.res \ $result/metaphlan2/lefse_Firmicutes.pdf # 批量绘制所有差异features柱状图 lefse-plot_features.py -f diff \ --archive none --format pdf \ temp/input.in temp/input.res \ $result/metaphlan2/lefse_
Kraken2可以快速完成读长(read)层面的物种注释和定量，还可以进行重叠群(contig)、基因(gene)、宏基因组组装基因组(MAG/bin)层面的序列物种注释。
# 方法1.启动kraken2工作环境 conda activate kraken2 # 方法2.启动指定位置的环境 # source /conda2/envs/kraken2/bin/activate # 记录软件版本 kraken2 --version # 2.1.1
{1}代表样本名字
输入：temp/qc/{1}_1_kneaddata_paired*.fastq 质控后的FASTQ数据
参考数据库：-db $/kraken2/mini/，默认标准数据库>50GB，这里使用8GB迷你数据库。
输出结果：每个样本单独输出，temp/kraken2/{1}_report和temp/kraken2/{1}_output
整合后的输出结果： result/kraken2/taxonomy_count.txt 物种丰度表
mkdir -p temp/kraken2
(可选) 单样本注释，5m
i=C1 # 1m，--use-mpa-style可直接输出metphlan格式，但bracken无法处理 # 2020/12/02版，65K双端序列，38.58%可注释，61.42%未注释，耗时5s，内存峰值8G # 2021/04/23版，65K双端序列，52.17%可注释，52.17%未注释，耗时5s，内存峰值8G # 内存8G的PC可运行，需要与硬盘交换，需3m，内存峰值4.5G kraken2 --db ${db}/kraken2/mini/ --paired temp/qc/${i}_?.fastq \ --threads 8 --use-names --report-zero-counts \ --report temp/kraken2/${i}.report \ --output temp/kraken2/${i}.output
多样本并行生成report，1样本8线程，内存大但速度快，内存不多不建议用多线程
tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \ "kraken2 --db ${db}/kraken2/mini --paired temp/qc/{1}_?.fastq \ --threads 8 --use-names --report-zero-counts \ --report temp/kraken2/{1}.report \ --output temp/kraken2/{1}.output"
使用krakentools转换report为mpa格式
for i in `tail -n+2 result/metadata.txt|cut -f1`;do kreport2mpa.py -r temp/kraken2/${i}.report \ --display-header \ -o temp/kraken2/${i}.mpa;done
合并样本为表格
mkdir -p result/kraken2 # 输出结果行数相同，但不一定顺序一致，要重新排序 tail -n+2 result/metadata.txt|cut -f1|rush -j 1 \ 'tail -n+2 temp/kraken2/{1}.mpa | LC_ALL=C sort | cut -f 2 | sed "1 s/^/{1}\n/" > temp/kraken2/{1}_count ' # 提取第一样本品行名为表行名 header=`tail -n 1 result/metadata.txt | cut -f 1` echo $header tail -n+2 temp/kraken2/${header}.mpa | LC_ALL=C sort | cut -f 1 | \ sed "1 s/^/Taxonomy\n/" > temp/kraken2/0header_count head -n3 temp/kraken2/0header_count # paste合并样本为表格 ls temp/kraken2/*count paste temp/kraken2/*count > result/kraken2/tax_count.mpa # 检查表格及统计 csvtk -t stat result/kraken2/tax_count.mpa
参数简介：
循环重新估计每个样品的丰度
# 设置估算的分类级别D,P,C,O,F,G,S，常用 P和S tax=P mkdir -p temp/bracken for i in `tail -n+2 result/metadata.txt|cut -f1`;do # i=C1 bracken -d ${db}/kraken2/mini \ -i temp/kraken2/${i}.report \ -r 100 -l ${tax} -t 0 \ -o temp/bracken/${i};done
结果描述：共7列，分别为物种名、ID、分类级、读长计数、补充读长计数、总数、百分比
name taxonomy_id taxonomy_lvl kraken_assigned_reads added_reads new_est_reads fraction_total_reads Capnocytophaga sputigena 1019 S 4798 996 5794 0.23041 Capnocytophaga sp. oral taxon 878 1316596 S 239 21 260 0.01034
bracken结果合并成表
# 输出结果行数相同，但不一定顺序一致，要去年表头重新排序 # 仅提取第6列reads count，并添加样本名 tail -n+2 result/metadata.txt|cut -f1|rush -j 1 \ 'tail -n+2 temp/bracken/{1} | LC_ALL=C sort | cut -f6 | sed "1 s/^/{1}\n/" > temp/bracken/{1}.count ' # 提取第一样本品行名为表行名 h=`tail -n1 result/metadata.txt|cut -f1` tail -n+2 temp/bracken/${h}|sort|cut -f1 | \ sed "1 s/^/Taxonomy\n/" > temp/bracken/0header.count # 检查文件数，为n+1 ls temp/bracken/*count | wc # paste合并样本为表格，并删除非零行 paste temp/bracken/*count > result/kraken2/bracken.${tax}.txt # 统计行列，默认去除表头 csvtk -t stat result/kraken2/bracken.${tax}.txt
结果筛选
# 需要指定安装R的位置和脚本位置 # alias Rscript="/anaconda2/bin/Rscript --vanilla" sd=/db/EasyMicrobiome/script # microbiome_helper按频率过滤，-r可标准化，-e过滤 Rscript $sd/filter_feature_table.R \ -i result/kraken2/bracken.${tax}.txt \ -p 0.01 \ -o result/kraken2/bracken.${tax}.0.01 # > 0.01(1%)的样本在出现，数量会明显减少 csvtk -t stat result/kraken2/bracken.${tax}.0.01 # 门水平去除脊索动物 grep 'Chordata' result/kraken2/bracken.P.0.01 grep -v 'Chordata' result/kraken2/bracken.P.0.01 > result/kraken2/bracken.P.0.01-H # 按物种名手动去除宿主污染，以人为例(需按种水平计算相关结果) # 种水平去除人类P:Chordata,S:Homo sapiens grep 'Homo sapiens' result/kraken2/bracken.S.0.01 grep -v 'Homo sapiens' result/kraken2/bracken.S.0.01 > result/kraken2/bracken.S.0.01-H
分析后清理每条序列的注释大文件
# rm -rf temp/kraken2/*.output
多样性分析/物种组成，详见3StatPlot.sh，Kraken2结果筛选序列见附录