【生信】使用QIIME进行 进化树，Alpha,Beta多样性 分析

使用QIIME进行 进化树，Alpha,Beta多样性 分析

---

上回讲到，使用Usearch进行进化树，Alpha,Beta多样性的分析。同时，我们还要再次强调QIIME的伟大之处在于全流程分析的能力。因此，使用QIIME同样可以进行上述分析过程。因此，将QIIME进行进化树，Alpha,Beta多样性 分析的脚本粘贴如下，怀着开放科研的心态，供有需要的朋友参考。同时，需要注意的是因为每个人的使用环境不同，可能部分无法正常运行，因此这对于新手来说可能不是那么友好。

#!/bin/bash
## 使用QIIME进行 进化树，Alpha,Beta多样性 分析

## 文件
HOME=`pwd`                ### 操作空间
temp=${HOME}/temp2              ### 临时文件存放目录
res=${HOME}/res2                ### 结果文件存放目录
data=${HOME}/seq                ### 测序数据存放位置
refdb=${HOME}/database          ### 参考测序数据库
err=${HOME}/err2                ### 错误信息
log=${HOME}/log2                ### 日志信息

printf "最终的OTU表:%s\t最终的OTU序列:%s\n"${res}/otu_table_tax.txt ${res}/rep_seqs.fasta
echo "进化树构建"
# clustalo多序列比对，如果没有请安装Clustal Omega
clustalo -i ${res}/rep_seqs.fasta -o ${temp}/rep_seqs_align.fa --seqtype=DNA --full --force --threads=30

echo "筛选结果中保守序列和保守区"
filter_alignment.py -i ${temp}/rep_seqs_align.fa -o ${temp}/  # rep_seqs_align_pfiltered.fa, only very short conserved region saved

echo "基于fasttree建树"
make_phylogeny.py -i ${temp}/rep_seqs_align_pfiltered.fasta -o ${res}/rep_seqs.tree # generate tree by FastTree

echo "Alpha多样性"
echo "查看样品的数据量最小值"
biom summarize-table -i ${res}/otu_table3.biom

echo "基于最小值进行重抽样标准化"
single_rarefaction.py -i ${res}/otu_table3.biom -o ${temp}/otu_table_rare.biom -d 2797

echo "计算常用的四种Alpha多样性指数"
alpha_diversity.py -i ${temp}/otu_table_rare.biom -o ${res}/alpha.txt -t ${res}/rep_seqs.tree -m shannon,chao1,observed_otus,PD_whole_tree

echo "Beta多样性"
echo "CSS标准化OTU表"
normalize_table.py -i ${res}/otu_table3.biom -o ${temp}/otu_table_css.biom -a CSS

echo "转换标准化OTU表为文本，用于后期绘图"
biom convert -i ${temp}/otu_table_css.biom -o ${res}/otu_table_css.txt --table-type="OTU table" --to-tsv

echo "删除表格多余信息，方便R读取"
sed -i '/# Const/d;s/#OTU //g;s/ID.//g' ${res}/otu_table_css.txt

echo "计算Beta多样性"
beta_diversity.py -i ${temp}/otu_table_css.biom -o ${res}/beta/ -t ${res}/rep_seqs.tree -m weighted_unifrac,unweighted_unifrac

echo "Beta多样性距离文件整理，方便R读取"
sed -i 's/^\t//g' ${res}/beta/*

echo "物种分类统计，筛选进化树和其他操作"
cd ${HOME}
echo "按物种分类级别分类汇总"
echo "结果按门、纲、目、科、属五个级别进行分类汇总，对应结果的L2-L6"
summarize_taxa.py -i ${res}/otu_table3.biom -o ${res}/sum_taxa # summary each level percentage
echo "修改一下文本表头，适合R读取的表格格式"
sed -i '/# Const/d;s/#OTU ID.//g' ${res}/sum_taxa/* # format for R read
#echo "以门为例查看结果"
#less -S ${res}/sum_taxa/otu_table3_L2.txt

echo "筛选可展示的进化树"
echo "选择OTU表中丰度大于0.1%的OTU"
filter_otus_from_otu_table.py --min_count_fraction 0.001 -i ${res}/otu_table3.biom -o ${temp}/otu_table_k1.biom

echo "获得对应的fasta序列"
filter_fasta.py -f ${res}/rep_seqs.fasta -o ${temp}/tax_rep_seqs.fa -b ${temp}/otu_table_k1.biom 

echo "统计序列数量，104条，一般100条左右即有大数据的B格，又能读懂和更清规律和细节"
grep -c '>' ${temp}/tax_rep_seqs.fa # 104

echo "多序列比对"
clustalo -i ${temp}/tax_rep_seqs.fa -o ${temp}/tax_rep_seqs_clus.fa --seqtype=DNA --full --force --threads=30

echo "建树"
make_phylogeny.py -i ${temp}/tax_rep_seqs_clus.fa -o ${temp}/tax_rep_seqs.tree

echo "格式转换为R ggtree可用的树"
sed "s/'//g" ${temp}/tax_rep_seqs.tree > ${res}/tax_rep_seqs.tree # remove '

echo "获得序列ID"
grep '>' ${temp}/tax_rep_seqs_clus.fa|sed 's/>//g' > ${temp}/tax_rep_seqs_clus.id

echo "获得这些序列的物种注释，用于树上着色显示不同分类信息"
awk 'BEGIN{OFS="\t";FS="\t"} NR==FNR {a[$1]=$0} NR>FNR {print a[$1]}' ${res}/rep_seqs_tax_assignments.txt ${temp}/tax_rep_seqs_clus.id|sed 's/; /\t/g'|cut -f 1-5 |sed 's/p__//g;s/c__//g;s/o__//g' > ${res}/tax_rep_seqs.tax

echo " 其他操作"
#echo "将mappingfile转换为R可读的实验设计"
#sed 's/#//' mappingfile.txt > ${res}/design.txt

echo "转换文本otu_table格式为R可读"
sed '/# Const/d;s/#OTU //g;s/ID.//g' ${res}/otu_table3.txt > ${res}/otu_table.txt

echo "转换物种注释信息为制表符分隔，方便R读取"
sed 's/;/\t/g;s/ //g' ${res}/rep_seqs_tax_assignments.txt > ${res}/rep_seqs_tax.txt