擴增子分析解讀2提取barcode,質控及樣品拆分,切除擴增引物
本網對Markdown排版支援較差，請跳轉“巨集基因組”公眾號閱讀；
寫在前面
之前釋出的《擴增子圖表解讀》系列，相信很多朋友都看過了(連結直達7月文章目錄)。
這些內容的初衷是寫給本領域剛進實驗室的學生讀，加速大家對同行文章的解讀能力。如果連同行的結果都看不懂，何談對資料的理解，對科學問題的解釋。希望剛入行的小夥伴多讀高水平文章，配合我的解讀，定能讓理解上升一個層次。
《擴增子分析解讀》系列文章介紹
擴增子分析是目前巨集基因組研究中最常用的技術，由於微生物組受環境影響大，實驗間重複較差，更需要更多的實驗重複和分析技術來保證結果的準確性、可重複性。
本系統文章叫分析解讀，即有詳細的分析流程程式碼，又有本人對使用引數、備選引數意義的解讀，可以讓大部分零基礎的人，更好的理解資料分析過程，並可親自實踐在自己的課題上，獲得更好、更合理的實驗結果。
16S擴增子測序資料主要來自HiSeq2500產出的雙端各250 bp (PE250)資料，因為讀長長且價格便宜(價效比高)。HiSeqX PE150和MiSeq PE300也比較常見，但PE150過短解析度低，而PE300價格高且末端序列質量過低。此外454在之前研究較但裝置已經停產，PacBio讀長長可直接測序16S全長1.5kb代表未來的趨勢。
本文采用目前最主流的擴增子測序資料型別HiSeq2500 PE250型別資料為例，結合目前主流方法QIIME+USearch定製的分析流程。本課程中所需的測序資料、實驗設計和課程分析生成的中間檔案，均可以直去百度雲下載。連結：http://pan.baidu.com/s/1hs1PXcw
本課程程式碼的執行，至少需要Linux平臺+安裝QIIME1.9.1，我之前釋出過三種安裝QIIME的方法詳見文章目錄，總有一款適合你。
第二節. 提取barcode,樣品拆分及質控,切除擴增引物
先看一下擴增子分析的整體流程，從下向上逐層分析。
分析前準備
# 進入工作目錄
cd example_PE250上一節，我們拿到了雙端資料，進行了質控、並對實驗設計進行了填寫和檢查無誤、最後將雙端資料合併為單個檔案進行下游分析。
接下來我們將序列末端的barcode標籤切下來，因為它們是人為添長的不屬於實驗物件；再根據標籤序列與實驗設計檔案比對，對每條序列屬於每個樣品進行分類；最後我們切除掉擴增使用的引物，因為每個序列均一樣，不切除反而容易引用錯誤。這樣我們就獲得了高質量的擴增區域資料，並且序列名稱中包括了樣品資訊。
4. 提取barcode
為什麼擴增子有barcode？
我以前做過sRNA-Seq, RNA-Seq, ChIP-Seq等等，都是一樣文庫一個樣品，為什麼沒有用過barcode。
因為擴增子目前研究物件細菌、真菌多樣性沒有表達基因數量大，一般是幾百到千的水平，對資料量要求最多10萬條序列即可飯飽合。而Illumina測序儀的通量很高，採用Index來區分每個文庫，仍然每個文庫的資料量達到千萬的級別，加上建庫測序的成本也不會低於千元。對於擴增子動輒成百上千的樣品即太貴，又浪費。因此將擴增子樣本新增上barcode(標籤)，通常將48/60個樣品混合在一起，構建一個測序文庫，達到高通量測序大量樣品同時降低實驗成本的目的。
通常的測序儀下機資料，只經過Index比對，拆分成來自不同文庫的資料檔案，分發給使用者。而擴增子的一個文庫包括幾十個樣品，還區要通過每個樣品上標記的特異Barcode進一步區分，再進行下游分析。
注：如果你是自己構建測序文庫，這步需要使用者手動拆分樣品；如果是公司建庫並測序，他們有樣品的barcode資訊，通常也會返給使用者樣品拆分後的資料，無需本節中的操作。但原理還是要理解，對整體分析的把握非常重要。
Barcode在擴增子的位置和型別？
Barcode位於引物的外側，比較典型的有三種，上圖展示的為最常用的barcode位於左端(上引物上游)，此外還有右端和雙端兩類也比較常用。
本分析採用的資料型別為右端barcode。
extract_barcodes.py是QIIME中用於切除barcode的指令碼，支援你所想到的所有型別。
-f引數為輸入檔案，即上文中合併雙端資料後的檔案；
-m為實驗設計檔案；
-o為輸出切下的barcode的資料目錄；
-c為barcode型別選擇，目前的barcode_paired_stitched選擇支援右端和雙端型別，如果是左端型別，請改為barcode_single_end；
–bc1_len設定左端barcode的長度，按實驗設計添寫，本資料只有右端則為零；
–bc2_len設定右端barcode的長度，按實驗設計添寫，本資料只有右端為6，即長度為6bp，常用的還有8，10；
-a是使用實驗設計中的引物來調整序列的方向，本實驗的測序無方向，必須開此選項調整，有些公司的建庫型別有方向，但開了總沒錯；
–rev_comp_bc2是將右端barcode取反向互補再與左端相連，建議開啟，這樣你的實驗設計中barcode一欄直接將添寫原始barcode即可，不用再考慮方向了；如果是雙端則將兩個barcode直接連起來填在barcode列，不用考慮方向。
# 提取barcode
extract_barcodes.py -f temp/PE250_join/fastqjoin.join.fastq \
-m mappingfile.txt \
-o temp/PE250_barcode \
-c barcode_paired_stitched --bc1_len 0 --bc2_len 6 -a --rev_comp_bc2這步任務會在輸出目錄temp/PE250_barcode中生成5個檔案
barcodes.fastq # 切下來的barcode，用於後序拆分樣品
barcodes_not_oriented.fastq # 方向不確定序列的barcode。連引物都不匹配，質量太差，建議不再使用
reads1_not_oriented.fastq # 方向不確定序列的序列，可能barcode切錯方向。連引物都不匹配，質量太差，建議不再使用
reads2_not_oriented.fastq # 空檔案
reads.fastq # 序列檔案，與barcode對應，用於下游分析5. 質控及樣品拆分
上步對序列方向進行了調整，並切開了barcode與擴增序列。下面使用split_libraries_fastq.py對混池根據barcode拆分樣品，同時篩選質量大於20(即準確度99%)的序列進行下游分析。引數解釋如下：
-i 輸入序列檔案；
-b 輸入barcode檔案，上步產生；
-m 實驗設計，依賴樣品barcode列表將序列重複命名
-q 測序質量閾值，20為99%準確率，30為99.9%準確
–min_per_read_length_fraction 最小高質量序列比例，0.75即只少有連序75%的鹼基質量高於20
–max_barcode_errors barcode允許的鹼基錯配個數，0為不允許，即使修改通常也不允許，不信你試試
barcode_type調置barcode型別，預設為golay_12，並支援錯配；我們通常添整數，為barcode的長度總和，本實驗中添0+6=6。
# 質控及樣品拆分
split_libraries_fastq.py -i temp/PE250_barcode/reads.fastq \
-b temp/PE250_barcode/barcodes.fastq \
-m mappingfile.txt \
-o temp/PE250_split/ \
-q 20 --max_bad_run_length 3 --min_per_read_length_fraction 0.75 --max_barcode_errors 0 --barcode_type 6
執行結果會輸出三個檔案
histograms.txt # 所有序列長度分佈資料
seqs.fna # 質控並拆分後的資料，序列按樣品編號為SampleID_0/1/2/3
split_library_log.txt # 日誌檔案，有基本統計資訊和每個樣品的資料量6. 切除引物序列
這裡我們介紹一款高通量引物切除軟體，cutadapt，2017-6-16最新版1.14；
https://pypi.python.org/pypi/cutadapt
此軟體發2011年釋出至今一直在更新，Google Scholar截止17年8月8日引用2263次。
Cutadapt 1.14下載和安裝
# 下載，請儘量檢查主頁下載最新版原始碼
wget https://pypi.python.org/packages/16/e3/06b45eea35359833e7c6fac824b604f1551c2fc7ba0f2bd318d8dd883eb9/cutadapt-1.14.tar.gz
# 解壓
tar xvzf cutadapt-1.14.tar.gz
# 進入程式目錄
cd cutadapt-1.14/
# 安裝在當前使用者目錄，不需管理員許可權
python setup.py install --usercutadapt切除雙端引物及長度控制，引數詳解：
-g 5’端引物
-a 3’端引物，需要將引物取反向互補
-e 引物匹配允許錯誤率，我調置0.15，一般引物20bp長允許3個錯配，為了儘量把引物切乾淨
-m 最小序列長度，根據情況設定，本實驗擴增V5-V7區，長度主要位於350-400，故去除長度小於300bp的序列，無經驗可不填此引數
–discard-untrimmed 如果引物末切掉的序列直接放棄
seqs.fna 為輸入檔案
PE250_P5.fa 為輸出檔案
cutadapt -g AACMGGATTAGATACCCKG -a GGAAGGTGGGGATGACGT -e 0.15 -m 300 --discard-untrimmed temp/PE250_split/seqs.fna -o temp/PE250_P5.fa程式執行結果會在螢幕上輸出結果統計報告，如去接頭比例、去掉過短序列比例等；還有去除引物的長度分佈，看看有沒有異常。如3’端序列沒有取反向互補會無法去除19bp序列，而是幾bp的錯誤序列。
Cutadapt結果報告
This is cutadapt 1.14 with Python 2.7.12
Command line parameters: -g AACMGGATTAGATACCCKG -a GGAAGGTGGGGATGACGT -e 0.15 -m 300 --discard-untrimmed temp/PE250_split/seqs.fna -o temp/PE250_P5.fa
Trimming 2 adapters with at most 15.0% errors in single-end mode ...
Finished in 41.03 s (32 us/read; 1.87 M reads/minute).
=== Summary ===
Total reads processed: 1,277,436
Reads with adapters: 1,277,194 (100.0%)
Reads that were too short: 8,849 (0.7%)
Reads written (passing filters): 1,268,345 (99.3%)
Total basepairs processed: 522,379,897 bp
Total written (filtered): 495,607,409 bp (94.9%)
=== Adapter 1 ===
Sequence: GGAAGGTGGGGATGACGT; Type: regular 3'; Length: 18; Trimmed: 202757 times.
No. of allowed errors:
0-5 bp: 0; 6-12 bp: 1; 13-18 bp: 2
Bases preceding removed adapters:
A: 96.3%
C: 1.5%
G: 0.8%
T: 1.3%
none/other: 0.0%
WARNING:
The adapter is preceded by "A" extremely often.
The provided adapter sequence may be incomplete.
To fix the problem, add "A" to the beginning of the adapter sequence.
Overview of removed sequences
length count expect max.err error counts
3 3 19959.9 0 3
4 4 4990.0 0 4
6 2 311.9 0 2
8 1 19.5 1 1
11 1 0.3 1 1
13 1 0.0 1 1
15 9 0.0 2 9
17 42 0.0 2 42
18 202626 0.0 2 202626
19 56 0.0 2 56
20 1 0.0 2 1
21 1 0.0 2 1
32 1 0.0 2 1
38 1 0.0 2 1
39 1 0.0 2 1
41 1 0.0 2 1
309 2 0.0 2 2
310 1 0.0 2 1
311 3 0.0 2 3
=== Adapter 2 ===
Sequence: AACMGGATTAGATACCCKG; Type: regular 5'; Length: 19; Trimmed: 1074437 times.
No. of allowed errors:
0-5 bp: 0; 6-12 bp: 1; 13-19 bp: 2
Overview of removed sequences
length count expect max.err error counts
3 2 19959.9 0 2
7 1 78.0 1 0 1
8 2 19.5 1 1 1
10 6 1.2 1 4 2
11 1 0.3 1 1
12 3 0.1 1 2 1
13 5 0.0 1 3 2
14 24 0.0 2 17 7
15 51 0.0 2 32 14 5
16 71 0.0 2 56 12 3
17 134 0.0 2 92 30 12
18 327 0.0 2 189 117 21
19 1059175 0.0 2 1056863 2069 243
20 13846 0.0 2 1817 10955 1074
21 744 0.0 2 5 10 729
22 1 0.0 2 1
23 2 0.0 2 2
24 1 0.0 2 1
25 2 0.0 2 2
27 5 0.0 2 5
28 2 0.0 2 2
29 2 0.0 2 2
30 1 0.0 2 1
31 2 0.0 2 2
32 10 0.0 2 10
49 1 0.0 2 1
51 1 0.0 2 1
166 1 0.0 2 1
291 6 0.0 2 6
401 2 0.0 2 2
409 1 0.0 2 1
443 1 0.0 2 1
460 2 0.0 2 2
465 2 0.0 2 2
WARNING:
One or more of your adapter sequences may be incomplete.
Please see the detailed output above.寫在後面
今天先到這裡，本文已經講了三個程式的使用，夠大家學習一會的了。要想了解這些程式的更多功能，一定要閱讀程式的幫助全文，才能有更深入的理解。
下節預告：擴增子分析解讀3去冗餘,聚類,生成OTU表
Reference
想了解更多16S/ITS/18S擴增子、巨集基因組、巨集轉錄組文獻閱讀和分析相關文章，快關注“巨集基因組”公眾號，乾貨第一時間推送。
系統學習生物資訊，快關注“生信寶典”，那裡有幾千志同道合的小夥伴一起學習。
