2.3 處理 發展 es2017 根據 條件 一個 命名 reads

本節課程，需要先完成 擴增子分析解讀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格式文件，用於後續分析

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