擴增子分析解讀4去嵌合體 非細菌序列 生成表明性序列和OTU表

本節課程，須要先完成

擴增子分析解讀1質控 實驗設計 雙端序列合併

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

3格式轉換 去冗餘 聚類

 

先看一下擴增子分析的總體流程，從下向上逐層分析

分析前準備

# 進入工做目錄
cd example_PE250

上一節回顧：咱們製做了Usearch要求格式的Fasta文件，對全部序列進行去冗餘和低丰度過濾，並聚類生成了OTU。

 

接下來咱們對OTU進一步去除嵌合體，並生成表明性序列和OTU表。

 

什麼是chimeras(嵌合體)？

嵌合體序列由來自兩條或者多條模板鏈的序列組成，示意圖以下：

在PCR反應中，延伸階段因爲不徹底延伸，就會致使嵌合體序列的出現，以上圖爲例，在擴增序列X的過程當中，在序列延伸階段，只產生了部分X序列延伸階段就結束了，在下一輪的PCR反應中，這部分序列做爲序列Y的引物接着延伸，擴增就會造成X和Y的嵌合體序列；

 

在放一張具體一點的示意圖，不徹底延伸產生的序列做爲下一輪PCR反應的產物，進行延伸

一般在PCR過程當中，大概有1%的概率會出現嵌合體序列，在16S/18S/ITS 擴增子測序的分析中，系統類似度極高，嵌合體可達1%-20%，須要去除嵌合體序列。

 

嵌合體的比例與PCR循環數相關，循環數越高，嵌合體比例越高。

 

有玩過魔獸有小夥伴記得精靈族的終極兵種雙頭龍奇美拉嗎？它的英文就是chimera，即中文的嵌合體，奇美拉是音譯。

10. 基於數據庫去嵌合體(可選)

上文第9步，聚OTU時，已經按照組內的序列類似狀況，直接denovo去除了大量嵌合體。目前這步基於數據庫去嵌合體，在之前的分析中是必作的，但隨着技術發展，發現這步可能也會形成假陰性。讀者能夠實驗設計、初步結果和預期來判斷是否須要這步處理。本文示例對每一步均進行操做，便是我的風格，又是爲了給你們展示一個比較全面的流程。以前Usearch做者推薦使用RDP數據去嵌合，並提供了下載連接；如今做者建議，若是作，就用Sliva或Unite這種全面的大數據庫，不推薦用RDP這種小數據庫，之前的建議是錯的。軟件方法均是不斷進步的，我尚未系統比較做者的新建議有多大改進，這裏仍是按照原來的方法進行，讀者能夠自行嘗試新方法。

# 下載Usearch推薦的參考數據庫RDP
wget http://drive5.com/uchime/rdp_gold.fa
# 基於RDP數據庫比對去除已知序列的嵌合體
./usearch10 -uchime2_ref temp/otus.fa \
 -db rdp_gold.fa \
 -chimeras temp/otus_chimeras.fa \
 -notmatched temp/otus_rdp.fa \
 -uchimeout temp/otus_rdp.uchime \
 -strand plus -mode sensitive -threads 96

採用-uchime2_ref參數去嵌合體，後面接OTU序列(輸入文件)；

-db 指定參考數據庫，這裏用RDP；

-chimeras 輸出檢測爲嵌合體的序列;

-notmatched 輸出不匹配數據庫的結果，即非嵌合，非相同序列；

-uchimeout 輸入嵌合體的檢測詳細信息，如每一個嵌合體的來源，與那幾個親本類似等；

-strand 指定鏈方向，通常爲正；

-mode 選擇模式，敏感的代價是嵌合體鑑定的高假陽性率；

-threads 設計線程數，程序默認系統小於10個線程爲單線程；多於10個線程爲10線程，根據實際狀況設置，不清楚不用更好。

 

上面計算結果Chimeras 2669/5489 (48.6%), in db 51 (0.9%), not matched 2769 (50.4%)，即5489個OTU有2669檢測爲嵌合、51個與數據庫序列一致爲非嵌合，另外2769與數據庫不匹配不肯定是否爲嵌合。對應temp/otus_rdp.uchime文件中第三列的Y/N/?

 

咱們想要是的排除嵌合的部分，即51+2769=2820。思路是將所有OTU中鑑定爲嵌合的排除掉。

# 得到嵌合體的序列ID
grep '>' temp/otus_chimeras.fa | sed 's/>//g' > temp/otus_chimeras.id
# 剔除嵌合體的序列
filter_fasta.py -f temp/otus.fa -o temp/otus_non_chimera.fa -s temp/otus_chimeras.id -n
# 檢查是否爲預期的序列數量2820
grep '>' -c temp/otus_non_chimera.fa

11. 去除非細菌序列(可選)

此步也是非必須的，容易形成假陰性。分析中有不少我的習慣的因素在裏面，因此不一樣人的分析結果，也會略有不一樣。也缺乏系統的評估到底那些更好，由於好與很差是有條件的，如何判斷也不容易説清楚，這就是經驗；項目經驗是通過大量的項目反覆研究積累出來的。

我的習慣在大數據面前，結果再多也沒用，得找到有意義的東西，因此原則上是能捨即舍，更容易發現規律。萬一沒有發現，再回去把扔掉的撿回來試試。若是什麼都不仍，規律可能永遠藏在大數據的海洋中。

 

這步的原理是將OTU與Greengene (http://greengenes.secondgenome.com)的Align數據庫比對，篩選序列類似性大於75%以上的序列做爲細菌序列；此步能夠排除外源非細菌的污染，非細菌序列在接下來的分析中沒法註釋物種分類，也很難分析。

# 下載Greengene最新數據庫，320MB
wget -c ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz
# 解壓數據包後大小3.4G
tar xvzf gg_13_8_otus.tar.gz
# 將OTU與97%類似聚類的表明性序列多序列比對，大約8min
time align_seqs.py -i temp/otus_non_chimera.fa -t gg_13_8_otus/rep_set_aligned/97_otus.fasta -o temp/aligned/
# 沒法比對細菌的數量
grep -c '>' temp/aligned/otus_non_chimera_failures.fasta # 1860
# 得到不像細菌的OTU ID
grep '>' temp/aligned/otus_non_chimera_failures.fasta|cut -f 1 -d ' '|sed 's/>//g' > temp/aligned/otus_non_chimera_failures.id
# 過濾非細菌序列
filter_fasta.py -f temp/otus_non_chimera.fa -o temp/otus_rdp_align.fa -s temp/aligned/otus_non_chimera_failures.id -n
# 看咱們如今還有多少OTU:975
grep '>' -c temp/otus_rdp_align.fa

通過這一步過濾，從2820非嵌合的OTU，只剩下975個與細菌類似的OTU，這種數量才更接近真相。有些研究常常搞幾千、幾萬的OTU，假陽性結果90%以上，你以爲意義何在，如何指導下游實驗。

 

對於真菌ITS/18S，通常不建議用Unite數據庫去嵌合，由於ITS/18S在全部真核生物中都有，有待物種註釋後進一步確認。

12. 產生表明性序列和OTU表

表明性序列(representative sequences)即爲肯定的最終版的OTU，相似於參考基因組/cDNA將爲索引的字典。而後將全部數據mapping於OTU上來肯定各物種的丰度。

 

OTU表，是每一個OTU在每樣品中的丰度值，本質上每種高通量測序結果，都會有一個相似的表，如RNA-Seq是基因表達與樣品的表

# 重命名OTU，這就是最終版的表明性序列，即Reference(可選，我的習慣)
awk 'BEGIN {n=1}; />/ {print ">OTU_" n; n++} !/>/ {print}' temp/otus_rdp_align.fa > result/rep_seqs.fa
# 生成OTU表
./usearch10 -usearch_global temp/seqs_usearch.fa -db result/rep_seqs.fa -otutabout temp/otu_table.txt -biomout temp/otu_table.biom -strand plus -id 0.97 -threads 10
# 結果信息 01:20 141Mb   100.0% Searching seqs_usearch.fa, 32.3% matched
# 默認10線程，用時1分20秒，有32.3%的序列匹配到OTU上；用30線程反而用時3分04秒，不是線程越多越快，分發任務也是很費時間的

如今咱們得到了OTU表，用less temp/otu_table.txt查看一下吧。同時還有biom可處理的標準json格式文件，用於後續分析