分析前準備 # 進入工作目錄 cd example_PE250 上一節回顧：我們獲得了OTU序列的進化分析、同時計算Alpha和Beta多樣性值。 本節是最後一節，我們對物種進行分類統計，篩選高豐度結果用於進化樹展示，和其它用於R統計分析的結果生成 19. 按物種分類級別分類彙總 OTU表中最重要的 ...
分析前準備
# 進入工作目錄 cd example_PE250上一節回顧：我們獲得了OTU序列的進化分析、同時計算Alpha和Beta多樣性值。 本節是最後一節，我們對物種進行分類統計，篩選高豐度結果用於進化樹展示，和其它用於R統計分析的結果生成 19. 按物種分類級別分類彙總 OTU表中最重要的註釋信息是物種註釋信息。通常的物種註釋信息分為7個級別：界、門、綱、目、科、屬、種。種是最小的級別，和OTU類似但有不相同。 我們除了可以比較樣品和組間OTU水平差異外，還可以研究不同類似級別上的差異，它們是否存在那些共同的變化規律。 按照註釋的級別進行分類彙總，無論是Excel還R操作起來，都是很麻煩的過程。這裡我們使用QIIME自帶的腳本summarize_taxa.py。
# 結果按門、綱、目、科、屬五個級別進行分類彙總，對應結果的L2-L6 summarize_taxa.py -i result/otu_table4.biom -o result/sum_taxa # summary each level percentage # 修改一下文本表頭，適合R讀取的表格格式 sed -i '/# Const/d;s/#OTU ID.//g' result/sum_taxa/* # format for R read # 以門為例查看結果 less -S result/sum_taxa/otu_table4_L2.tx以門為例，我們看到樣品的OTU分佈在19個門，及每個門在各樣品中的相對比例。其它的各級別，用戶自己看吧。 這步的結果將用於後期統計和繪圖。 20. 篩選可展示的進化樹 我們在文章中看到幾種漂亮的進化樹，但是OTU通常成百上千，如果直接展示是根本看不清也是極醜的。 下麵教大家一些通常的方法來篩選數據，用於生成漂亮的進化樹。
# 選擇OTU表中豐度大於0.1%的OTU
filter_otus_from_otu_table.py --min_count_fraction 0.001 -i result/otu_table4.biom -o temp/otu_table_k1.biom
# 獲得對應的fasta序列
filter_fasta.py -f result/rep_seqs.fa -o temp/tax_rep_seqs.fa -b temp/otu_table_k1.biom
# 統計序列數量，104條，一般100條左右即有大數據的B格，又能讀懂和更清規律和細節
grep -c '>' temp/tax_rep_seqs.fa # 104
# 多序列比對
clustalo -i temp/tax_rep_seqs.fa -o temp/tax_rep_seqs_clus.fa --seqtype=DNA --full --force --threads=30
# 建樹
make_phylogeny.py -i temp/tax_rep_seqs_clus.fa -o temp/tax_rep_seqs.tree
# 格式轉換為R ggtree可用的樹
sed "s/'//g" temp/tax_rep_seqs.tree > result/tax_rep_seqs.tree # remove '
# 獲得序列ID
grep '>' temp/tax_rep_seqs_clus.fa|sed 's/>//g' > temp/tax_rep_seqs_clus.id
# 獲得這些序列的物種註釋，用於樹上著色顯示不同分類信息
awk 'BEGIN{OFS="\t";FS="\t"} NR==FNR {a[$1]=$0} NR>FNR {print a[$1]}' result/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' > result/tax_rep_seqs.tax
21. 其它
其它都是一些簡單的格式轉換，為後面統計分析而準備文件。
# 將mappingfile轉換為R可讀的實驗設計 sed 's/#//' mappingfile.txt > result/design.txt # 轉換文本otu_table格式為R可讀 sed '/# Const/d;s/#OTU //g;s/ID.//g' result/otu_table4.txt > result/otu_table.txt # 轉換物種註釋信息為製表符分隔，方便R讀取 sed 's/;/\t/g;s/ //g' result/rep_seqs_tax_assignments.txt > result/rep_seqs_tax.txt
