擴增子分析解讀5物種註釋 OTU表操做

本節課程，須要先完成《擴增子分析解讀》系列以前的操做

1質控 實驗設計 雙端序列合併

2提取barcode 質控及樣品拆分 切除擴增引物

3格式轉換 去冗餘 聚類

4去嵌合體 非細菌序列 生成表明性序列和OTU表

分析前準備

# 進入工做目錄
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.0

OE3: 2353.0

OE8: 3091.0

OE2: 3173.0

OE1: 3337.0

OE5: 3733.0

OE6: 4289.0

OE9: 4648.0

OE7: 5185.0

WT3: 10741.0

WT8: 12117.0

WT6: 14316.0

WT2: 14798.0

WT7: 14851.0

KO1: 14926.0

WT9: 15201.0

WT1: 15422.0

WT5: 15773.0

WT4: 16708.0

KO2: 17607.0

KO6: 23949.0

KO5: 26570.0

KO8: 27250.0

KO4: 32303.0

KO7: 33086.0

KO9: 35913.0

KO3: 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