本節課程，需要先完成《擴增子分析解讀》系列之前的操作
分析前准備
# 進入工作目錄 cd example_PE250
上一節回顧：我們學習了嵌合體的形成，以及基於參考數據庫去嵌合體；也學習了基於數據庫比對來篩選細菌或真菌；最后基於最確定的OTU，我們生成代表性序列和OTU表，這是每種高通量測序都有的結果，后續的結果將全部基於這兩個文件。
接下來我們學習對OTU進行物種注釋；OTU的操作，包括格式轉換、篩選添加物種信息、數據量篩選樣品、篩選高豐度的OTU、物種篩選等。
OTU表常用的BIOM格式
主頁：http://biom-format.org/ 。BIOM是英文The Biological Observation Matrix的縮寫，中文翻譯為生物觀測矩陣，是一種通過格式，用於生物學樣品對應觀測值的表格。它主要采用json/HD5F文件格式標准，即多維散列結構，保存表格結構數據結果。目前主流的宏基因組軟件均支持此格式文件，如QIIME、MG-RAST、PICRUSt、Mothur、phyloseq、MEGAN、VAMPS、metagenomeSeq、Phinch、RDP Classifier、USEARCH、PhyloToAST、EBI Metagenomics、GCModeller、MetaPhlAn 2。知道它有多重要了吧。
Biom文件處理系統biom程序是QIIME的必裝包，如果沒有安裝好，可嘗試下面步驟重裝
# 安裝依賴包 pip install numpy # 安裝biom格式轉換包 pip install biom-format # 安裝2.0格式支持 pip install h5py # 測序程序是否安裝成功 biom
13. 物種注釋
對於擴增子分析，最重要的就是物種信息。我們基於上節分析得到的代表性序列，采用上次已經下載的greengene的參考序列和物種注釋信息，比對軟件選擇rdp方法，進行注釋。
# 物種注釋 assign_taxonomy.py -i result/rep_seqs.fa \ -r gg_13_8_otus/rep_set/97_otus.fasta \ -t gg_13_8_otus/taxonomy/97_otu_taxonomy.txt \ -m rdp -o result
注：如果是ITS/18S數據，建議數據庫更改為UNITE，方法改為blast。詳細使用說明，請讀官方文檔http://qiime.org/scripts/assign_taxonomy.html
14. OTU表統計、格式轉換、添加信息
將OTU表轉換為Biom格式，這樣便於其它軟件對其操作。可添加上面獲得的物種信息，這樣表格的信息就更豐富了，再轉換為文本，便於人類可讀，同時使用summarize-table查看OTU表的基本信息。
# 文本OTU表轉換為BIOM：方便操作 biom convert -i temp/otu_table.txt \ -o result/otu_table.biom \ --table-type="OTU table" --to-json # 添加物種信息至OTU表最后一列，命名為taxonomy biom add-metadata -i result/otu_table.biom \ --observation-metadata-fp result/rep_seqs_tax_assignments.txt \ -o result/otu_table_tax.biom \ --sc-separated taxonomy --observation-header OTUID,taxonomy # 轉換biom為txt格式，帶有物種注釋：人類可讀 biom convert -i result/otu_table_tax.biom -o result/otu_table_tax.txt --to-tsv --header-key taxonomy # 查看OTU表的基本信息：樣品，OUT數量統計 biom summarize-table -i result/otu_table_tax.biom -o result/otu_table_tax.sum
現在我們獲得了OTU表的基本統計信息，用less result/otu_table_tax.sum查看一下吧，內容如下：
Num samples: 27 # 樣品數據Num observations: 975 # OTU數據Total count: 409647 # 總數據量Table density (fraction of non-zero values): 0.464 # 非零的單元格Counts/sample summary:Min: 2352.0 # 樣品數據量最小值Max: 35955.0 # 樣品數據量最大值Median: 14851.000 # 樣品數據量中位數Mean: 15172.111 # 樣品數據量平均數Std. dev.: 10691.823 # 樣品數據量標准變異Sample Metadata Categories: None provided # 樣品分類信息：末提供Observation Metadata Categories: taxonomy # 觀察值分類：物種信息Counts/sample detail: # 每個樣品的數據量OE4: 2352.0OE3: 2353.0OE8: 3091.0OE2: 3173.0OE1: 3337.0OE5: 3733.0OE6: 4289.0OE9: 4648.0OE7: 5185.0WT3: 10741.0WT8: 12117.0WT6: 14316.0WT2: 14798.0WT7: 14851.0KO1: 14926.0WT9: 15201.0WT1: 15422.0WT5: 15773.0WT4: 16708.0KO2: 17607.0KO6: 23949.0KO5: 26570.0KO8: 27250.0KO4: 32303.0KO7: 33086.0KO9: 35913.0KO3: 35955.0
biom的詳細使用說明，可以biom查看具體的功能，如添加注釋功能biom add-metadata --help可查看詳細說明。也可閱讀官網http://biom-format.org/
15. OTU表篩選
實驗中會有各種影響因素，我們要綜合各種背景知識來判斷如何篩選數據表，起到去偽存真，去粗取粗，由此及彼，有表及理的來回答科學問題。數據篩選是會運行分析流程和數據分析師的分水嶺。
看上面的的統計結果，樣本數據量從2k-35k，我們應去除過小的數據量樣本，提供更可能高的樣品最低豐度的數據用於下游標准化分析。這里我們選擇只保留數據量大於3000的樣品。
# 按樣品數據量過濾：選擇counts>3000的樣品 filter_samples_from_otu_table.py -i result/otu_table_tax.biom -o result/otu_table2.biom -n 3000 # 查看過濾后結果：只有25個樣品，975個OTU biom summarize-table -i result/otu_table2.biom
同時還要過濾低豐度的OTU，一般低於萬分之一豐度的菌，在功能研究可能還是比較困難的(早期文章454測序數據量少，通常只關注豐度千分之五以上的OTU)。
# 按OTU豐度過濾：選擇相對豐度均值大於萬分之一的OTU filter_otus_from_otu_table.py --min_count_fraction 0.0001 -i result/otu_table2.biom -o result/otu_table3.biom # 查看過濾后結果：只有25個樣品，346個OTU biom summarize-table -i result/otu_table3.biom
有些研究手段在特定有實驗中存在偏差，如2012Nature報導V5-V7在植物中擴增會偏好擴增Chloroflexi菌門，建議去除。
# 按物種篩選OTU表：去除p__Chloroflexi菌門 filter_taxa_from_otu_table.py -i result/otu_table3.biom -o result/otu_table4.biom -n p__Chloroflexi # 查看過濾后結果：只有25個樣品，307個OTU biom summarize-table -i result/otu_table4.biom
以上過濾條件是根據經驗、相關文獻設計的，如果不清楚，也不要隨便過濾，容易引起假陰性。
得到的最終結果，還要轉換為文本格式，和提取OTU表對應的序列，用於下游分析。
# 轉換最終biom格式OTU表為文本OTU表格 biom convert -i result/otu_table4.biom -o result/otu_table4.txt --table-type="OTU table" --to-tsv # OTU表格式調整方便R讀取 sed -i '/# Const/d;s/#OTU //g;s/ID.//g' result/otu_table4.txt # 篩選最終OTU表中對應的OTU序列 filter_fasta.py -f result/rep_seqs.fa -b result/otu_table4.biom -o result/rep_seqs4.fa