鉴定物种丰度

MetaPhlAn4

MetaPhlAn 4.0 tutorial

MetaPhlAn4是一个从宏基因组shotgun测序数据中分析微生物群落组成的工具

安装metaphlan

数据库位置 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

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_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% ，以此类推。

查看所有目级别的分类总和是100%

grep o__ SRS014476-Supragingival_plaque_profile.txt | grep -v f__ | cut -f1,3

使用toy数据库测试metaphlan

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.

2.1 准备HUMAnN2输入文件

小技巧：循环批量处理样本列表

    # 基于样本元数据提取样本列表命令解析
    # 去掉表头
    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减少输出文件加速

2.2 HUMAnN2计算物种和功能组成

• 物种组成调用MetaPhlAn2, bowtie2比对至核酸序列，解决有哪些微生物存在的问题；
• 功能组成为humann2调用diamond比对至蛋白库11Gb，解决这些微生物参与哪些功能通路的问题；
• 输入文件：temp/concat/*.fq 每个样品质控后双端合并后的fastq序列
• 输出文件：temp/humann2/ 目录下
  • C1_pathabundance.tsv
  • C1_pathcoverage.tsv
  • C1_genefamilies.tsv
• 整合后的输出：
  • result/metaphlan2/taxonomy.tsv 物种丰度表
  • result/metaphlan2/taxonomy.spf 物种丰度表（用于stamp分析）
  • result/humann2/pathabundance_relab_unstratified.tsv 通路丰度表
  • result/humann2/pathabundance_relab_stratified.tsv 通路物种组成丰度表
  • stratified(每个菌对此功能通路组成的贡献)和unstratified(功能组成)

启动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

2.3 物种组成表

2.3.1 样品结果合并

    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

2.3.2 转换为stamp的spf格式

    metaphlan_to_stamp.pl result/metaphlan2/taxonomy.tsv \
      > result/metaphlan2/taxonomy.spf
    head -n5 result/metaphlan2/taxonomy.spf
    # 下载metadata.txt和taxonomy.spf使用stamp分析
    # 网络分析见附录 metaphlan2-共有或特有物种网络图

2.3.3 (可选)Python绘制热图

    # 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

2.4 功能组成分析

2.4.1 功能组成合并、标准化和分层

合并通路丰度(pathabundance)，含功能和对应物种组成。
可选基因家族(genefamilies 太多)，通路覆盖度(pathcoverage)。
注意看屏幕输出# 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进行统计分析

2.4.2 差异比较和柱状图

两样本无法组间比较，在pcl层面替换为HMP数据进行统计和可视化。

参考 https://bitbucket.org/biobakery/humann2/wiki/Home#markdown-header-standard-workflow

• 输入数据：通路丰度表格 result/humann2/pathabundance.tsv
• 输入数据：实验设计信息 result/metadata.txt
• 中间数据：包含分组信息的通路丰度表格文件 result/humann2/pathabundance.pcl
• 输出结果：result/humann2/associate.txt

在通路丰度中添加分组

    ## 提取样品列表
    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

2.4.3 转换为KEGG注释

需要下载utility_mapping数据库并配置成功才可以使用。详见软件和数据库安装1soft_db.sh。

支持GO、PFAM、eggNOG、level4ec、KEGG的D级KO等注释，详见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等操作

2.5 GraPhlAn图

    # 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.6 LEfSe差异分析物种

• 输入文件：物种丰度表result/metaphlan2/taxonomy.tsv
• 输入文件：样品分组信息 result/metadata.txt
• 中间文件：整合后用于LefSe分析的文件 result/metaphlan2/lefse.txt，这个文件可以提供给www.ehbio.com/ImageGP 用于在线LefSE分析
• LefSe结果输出：result/metaphlan2/目录下lefse开头和feature开头的文件

前面演示数据仅有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_

2.7 Kraken2物种注释

Kraken2可以快速完成读长(read)层面的物种注释和定量，还可以进行重叠群(contig)、基因(gene)、宏基因组组装基因组(MAG/bin)层面的序列物种注释。

    # 方法1.启动kraken2工作环境
    conda activate kraken2
    # 方法2.启动指定位置的环境
    # source /conda2/envs/kraken2/bin/activate
    # 记录软件版本
    kraken2 --version # 2.1.1

2.7.1 Kraken2物种注释

• {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

2.7.2 Bracken估计丰度

参数简介：

• -d为数据库，与kraken2一致
• -i为输入kraken2报告文件
• r是读长，此处为100，通常为150
• l为分类级，本次种级别(S)丰度估计，可选域、门、纲、目、科、属、种：D,P,C,O,F,G,S
• t是阈值，默认为0，越大越可靠，但可用数据越少
• -o 输出重新估计的值

循环重新估计每个样品的丰度

    # 设置估算的分类级别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结果筛选序列见附录