擴增子分析解讀3格式轉換 去冗餘 聚類

本節課程，須要完成擴增子分析解讀1質控 實驗設計 雙端序列合併和2提取barcode 質控及樣品拆分 切除擴增引物

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

分析前準備

# 進入工做目錄
cd example_PE250

上一節回顧：咱們提取barcode，質控及樣品拆分，切除擴增引物，經歷了兩節課6步數據處理纔拿到咱們擴增的高質量目的片斷(貌似基因組/RNA-Seq測序結果直接就是這個階段了，能夠直接mapping)

 

接下來咱們將這些序列去冗餘、聚類爲OTU、再去除嵌合體，這樣就能夠得到高質量的OTU(相似於參考基因組/轉錄組)，用於定量分析每一個OTU的丰度。這一階段咱們使用著名的擴增子分析流程Usearch。

 

Usearch簡介

http://drive5.com/usearch/

Usearch以前介紹過http://www.cnblogs.com/freescience/p/7270861.html

軟件做者不只有Usearch一款軟件，它的Muscle(多序列比對，引用18659+4212次)，Uparse(OTU聚類算法，引用1529次), Uchime(擴增子嵌合體檢測，引用3558次)等衆多流行工具，我的引用超4萬次，並且發的軟件大多由做者一人完成，佩服。

 

Usearch安裝

這個軟件64位版收費，但32位對任何人免費，能夠在下面網址下載 http://www.drive5.com/usearch/download.html 贊成許可協議，選擇軟件版本(5.2 — 10.0)，選擇運行平臺(Linux, Windows或Mac OSX)填寫郵箱得到下載地址。不容許私人傳播。

這裏我選擇10.0版本，系統選擇Linux。

收到的郵件中第一個連接即下載地址，後面兩個連接爲幫助文檔和安裝說明，先不用管，按我下面的操做來。

# 下載程序並重命名：下載連接來自郵件，請用戶自行復制郵件中地址替換下面代碼中的網址；或者在windows裏面下載並重命名爲usearch10
wget -O "usearch10" http://drive5.com/cgi-bin/upload3.py?license=XXXXXX
# 添加可執行權限
chmod +x usearch10
# 運行程序測試，成功可顯示程序版本、系統信息和用戶受權信息
./usearch10

7. 格式轉換

作生信爲何要學Python/Perl/Shell這些語言，主要緣由是各軟件間要求的具體格式不一樣，須要進行格式轉換，才能繼續運行。所以想成爲高手，不會語言基本步履維艱。

 

咱們如今將QIIME拆分的結果類型，要轉換成Usearch要求的格式。常見的解決思路是讀Usearch幫助看它的格式要求，寫個Python/Perl腳本轉換格式。我這裏使用了Shell腳本一行解決，優勢是快，但缺點不少(人不容易看懂、不一樣Linux系統shell版本不一樣可能失效)

 

咱們要轉換的序列文件其實一直是fasta格式，只是序列名稱行格式不一樣

# 目前格式

>KO1_0 HISEQ:419:H55JGBCXY:1:1101:1931:2086 1:N:0:CACGAT orig_bc=TAGCTT new_bc=TAGCTT bc_diffs=0   

# Usearch要求的格式

>KO1_0;barcodelabel=KO1;

# 格式轉換
sed 's/ .*/;/g;s/>.*/&&/g;s/;>/;barcodelabel=/g;s/_[0-9]*;$/;/g' temp/PE250_P5.fa > temp/seqs_usearch.fa

上面這條命令有點複雜。sed是linux的一條命令，又是一種語言，擅長文本替換。替換的思路分四步：首先s/ ./;/g將原文件空格後面的內容(全是無用信息)替換爲分號；其次s/>./&&/g是將序列名重複一次；再次s/;>/;barcodelabel=/g將重複後的;>替換爲;barcodelabel=；最後s/_[0-9]*;$/;/g替換序列編號爲分號。這只是個人思路，分析數據如解答數學題，能夠有多種解法，你夠聰明還會想出更好的解法。

新人必定感受這命令每句都不像人話，我告訴你Perl和Shell就是這樣—難讀但高效。改用易讀的Python語言，確定沒有Shell簡潔。

 

8. 去冗餘

爲何要去冗餘？

由於原始序列幾百萬條，聚類計算的時間極其恐怖。而已知擴增子測序結果中序列重複度高，而且大量出現1次或幾回的序列統計學和功能上意義不大。所以將幾百萬條序列去冗餘，並過濾低丰度序列，通常只剩幾萬條，極大的減小了下游分析的工做量，並可以使結果更容易理解。

usearch10的去冗餘命令叫-fastx_uniques，緊跟着輸入文件；

-fastaout 接輸出文件；

-minuniquesize 參數設置保留的最小丰度reads數，建議最小設置爲2，去掉全部的單次出現序列(singletons)，數據量大建議設置總數據量的百萬分之一並取整數部分

-sizeout 在序列名稱中添加序列出現的頻率

# 序列去冗餘
./usearch10 -fastx_uniques temp/seqs_usearch.fa -fastaout temp/seqs_unique.fa -minuniquesize 2 -sizeout

計算過程當中出現以下信息：

00:06 607Mb   100.0% Reading temp/seqs_usearch.fa

00:06 574Mb  CPU has 96 cores, defaulting to 10 threads

00:08 915Mb   100.0% DF

00:09 935Mb  1268345 seqs, 686530 uniques, 624363 singletons (90.9%)

00:09 935Mb  Min size 1, median 1, max 18774, avg 1.85

62167 uniques written, 182874 clusters size < 2 discarded (26.6%)

主要內容爲讀取輸入文件；

檢查到系統有96個CPU，默認使用了10個線程；

總共有1268345條序列，其中非重複的序列有686530個，非重複且只出現一次的有624363個(90.9%的非冗餘序列是singletons，多嗎？)；

最小值、中位數、最大值、平均值；輸出結果有62167個結果，丟棄掉的數據佔26.6%。

 

本條命令的詳細使用，請閱讀官方文檔 http://www.drive5.com/usearch/manual/cmd_fastx_uniques.html

 

9. 聚類OTU

 

爲何要聚類OTU？

是由於Unique的序列仍然遠多於物種數量，而且擴增的物種可能存在rDNA的多拷貝且存在變異而獲得來自同一物種的多條序列擴增結果。目前人爲定義序列類似度一般97%以上爲OTU，大約是物種分類學種的水平，實際上1個OTU可能包括多個物種，而一個物種也可能擴增出多個OTU。

 

下面咱們用usearch10將非冗餘的序列聚類

-cluster_otus接輸入文件；

-otus後面爲輸出的otu文件的fasta格式；

-uparseout輸出聚類的具體細節

-relabel Otu爲重命名序列以Otu起始

# 聚類OTU
./usearch10 -cluster_otus temp/seqs_unique.fa -otus temp/otus.fa -uparseout temp/uparse.txt -relabel Otu

程序運行過程會顯示運行時間、進度，發現的OTU，以及嵌合體數據；結果以下：

04:11 84Mb    100.0% 5489 OTUs, 9209 chimeras

程序一共運行了3分39秒，聚類發現5486個OTUs，同時發現了9187個嵌合體並已被丟棄。

Usearch聚類算法之因此能發表在Nature Method上，就是由於其算法UParse在很是強的嵌合體檢測能力，對人工重組數據評估，更接近真實結果。下一節咱們將詳細講嵌合體產生的緣由，以及去除的原理。

 

本條命令的詳細使用，請閱讀官方文檔 http://www.drive5.com/usearch/manual/cmd_cluster_otus.html

 

小技巧：統計fasta文件中序列的數量

fasta文件每條序列以大於號(>)開始，其數量與序列數量相同，使用grep檢索含有>的行，同時用-c參數對數量進行統計，便可快速得到fasta文件中序列數量。

# 查看OTU數量
grep '>' -c temp/otus.fa