Metagenome多样性分析&差异分析 2023.01.05-01.07

多样性分析&差异分析

## 准备输入biom文件
ln -s ../01.taxon/out_S/S.biom
## 进行alpha多样性指数计算及绘制稀释曲线
singularity exec ../../software/MetaGenome Rscript ../script/diversity_alpha.R S.biom S.alpha
**output：
S.alpha.alpha-diversity.table alpha 多样性指数
S.alpha.Rarefaction_ggplot2.pdf 稀释曲线
S.alpha.Rarefaction_orginal.pdf 稀释曲线

## 进行beta多样性分析
singularity exec ../../software/MetaGenome Rscript ../script/diversity_beta.R S.biom sample.txt S.beta
**output：
S.beta.dist-braycurtis.table 距离矩阵
S.beta.NMDS.pdf NMDS-plot
S.beta.NMDS.table NMDS表格
S.beta.PCA.pdf PCA-plot
S.beta.PCA.table PCA表格
S.beta.PCoA.pdf PCoA-plot
S.beta.PCoA.table PCoA表格

差异分析

## 准备输入文件
ln -s ../01.taxon/sample.txt
ln -s ../01.taxon/out_S/S.taxName.count.tsv
## 数据准备-data部分
awk -F "\t" -v 'OFS=\t' '{$1=$NF; $NF=""; {print $0}}' \
S.taxName.count.tsv | \
sed 1d | \
sed 's/; /|/g' > S.lefse.tmp
## 数据准备-class部分
head -n 1 S.lefse.tmp |sed 's/\t/\n/g' | sed '1d' | \
while read sp ;do awk '$2=="'$sp'" {printf "\t"$1 }' sample.txt ;done |\
awk '{print "class"$0 }' > S.lefse.header
## 合并
cat S.lefse.header S.lefse.tmp > S.lefse
## 格式转换
singularity exec ../../software/lefse_latest.sif format_input.py \
S.lefse \ # 输入
S.lefse.in \ # 输出
-c 1 \ # class所在行
-s -1 \ # subclass所在行
-u 2 \ # subject所在行
-o 1000000 # normalization 到1M
## 运行lefse
singularity exec ../../software/lefse_latest.sif run_lefse.py \
S.lefse.in \ # 输入
S.lefse.out \ # 输出
-l 2 # LDA阈值
## LDA图
singularity exec ../../software/lefse_latest.sif plot_res.py \
S.lefse.out \
S.lefse.LDA.pdf \
--format pdf
## 进化分支图
singularity exec ../../software/lefse_latest.sif plot_cladogram.py \
S.lefse.out \
S.lefse.cladogram.pdf \
--format pdf \
--labeled_start_lev 1

也可以用DESeq2/edgeR 这类差异分析软件进行差异物种分析

## 准备count表格
ln -s ../01.taxon/out_S/S.count.tsv
## 替换掉空格
sed 's/ /_/g' S.count.tsv > S.count.matrix
## 进行差异分析
singularity exec ../../software/MetaGenome.sif perl \
../script/run_DE_analysis.pl \
--matrix S.count.matrix \ # count表格
--method DESeq2 \ # 分析方法，可设置为edgeR
--samples_file sample.txt \ # 样品分组信息
--output DESeq2_out # 输出结果目录
**output：
S.count.matrix.A_vs_B.DESeq2.DE_results ; 可基于FDR 和 Foldchange 进行差异物种筛选。常用标准为：
FDR < 0.05 && |log2(FC)|>1