基于qiime2 扩增子分析结果进行LEfse差异物种进化关系和出图

LEfSe（LDA Effect Size）分析是一种将非参数的Kruskal-Wallis以及Wilcoxon秩和检验，与线性判别分析（Linear discriminant analysis，LDA）效应量（Effect size）相结合的分析手段。与MetagenomeSeq分析类似，LEfSe分析也是一种差异分析方法；但LEfSe分析可以直接对门/纲/目/科/属/种的各级分类水平同时进行统计检验和差异分析。同时，LEfSe更强调寻找分组之间稳健的差异物种，即标志物种（Biomarker）。它的一大特点是，不仅局限于对不同样本分组中的群落组成差异进行分析，更可以深入到不同的子分组（Subgroup）中，挑取在不同子分组中表现一致的标志微生物类群。

要基于 QIIME 2 的扩增子分析结果进行 Lefse 分析，你需要将 QIIME 2 生成的结果文件转换成 Lefse 能够处理的格式，然后使用 Lefse 工具进行统计分析。以下是大致的步骤：

步骤概述：

1. 导出 QIIME 2 结果文件
2. 转换结果文件格式为 Lefse 可处理格式
3. 运行 Lefse 统计分析

环境安装

# bioconda安装
mamba create -n lefse -c bioconda lefse -y 

# 激活lefse环境
conda activate lefse 

# 查看输入数据格式化脚本
which lefse-format_input.py

转化最终结果到Level6

qiime taxa collapse \
    --i-table rst_table.qza \
    --o-collapsed-table collapse.table.qza \
    --p-level 6 \
    --i-taxonomy final_tax_sliva.qza

统计计算相对丰度

qiime feature-table relative-frequency \
    --i-table collapse.table.qza \
    --o-relative-frequency-table collapse.frequency.table.qza \
    --output-dir ./ 

导出biom文件

qiime tools export \
    --input-path collapse.frequency.table.qza \
    --output-path ./

将biom格式转换为普通文本

biom convert \
    -i feature-table.biom \
    -o collapse.frequency.table.tsv \
    --header-key "taxonomy" \
    --to-tsv

按tax水平过滤物种

sed 's/;/\|/g' collapse.frequency.table.tsv | \
    awk '{split($1, a, "|");if( a[6] != "__"){print $0}}' | \
    #sed 's/d\_\_Bacteria|//g' | \
    grep -vE "g__uncultured|d__Archaea|p__WPS-2|p__SAR324_clade|Constructed" | \
    sed 's/#OTU ID/Group/g;s/taxonomy//g' > collapse.frequency.table.lefse.tsv

运行lefse分析

lefse分析的输入文件大概格式：

前面两行是处理和分组，第三行是样品编号，第四行开始就是物种丰度了，第一列是物种名当然也可以是其他的，只要格式一致的数据都可以使用这个差异分析。

mamba activate lefse
# convert text file into lefse.input file 
# 数据输入format
lefse_format_input.py \
    collapse.frequency.table.lefse.tsv \
    collapse.frequency.table.lefse.in \
    -c 1 \
    -m f \
    -o 100000

# run lefse
run_lefse.py \
    collapse.frequency.table.lefse.in \
    collapse.frequency.table.lefse.res 

# select significant result Lefse
grep -E "HTN|Normal" \
        collapse.frequency.table.lefse.res \
        > collapse.frequency.table.lefse_signif.res

# plot lda 
lefse_plot_res.py \
    collapse.frequency.table.lefse_signif.res \
    lefse_final_lda.pdf \
    --format pdf \
    --autoscale 0

# plot cladogram 
lefse_plot_cladogram.py \
    collapse.frequency.table.lefse_signif.res \
    lefse_total_clado.pdf \
    --format pdf