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

本節課程，須要完成 擴增子分析解讀1質控 實驗設計 雙端序列合併

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

分析前準備

# 進入工做目錄
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的長度，按實驗設計添寫，本數據右端長度爲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對應，用於下游分析

更多說明建議閱讀幫助 http://qiime.org/scripts/extract_barcodes.html

 

5. 質控及樣品拆分

上步對序列方向進行了調整所有爲正向，並切開了barcode與擴增序列。下面使用split_libraries_fastq.py對混池根據barcode拆分樣品，同時篩選質量大於20(即準確度99%)的序列進行下游分析。參數解釋以下：

-i 輸入序列文件，上步產生；

-b 輸入barcode文件，上步產生；

-m 實驗設計，依賴樣品barcode列表將序列按樣品從新命名

-q 測序質量閾值，20爲99%準確率，30爲99.9%準確

--max_bad_run_length 容許連續的低質量鹼基調用的最大值，默認值爲3

--min_per_read_length_fraction 最小高質量序列比例，0.75即最少有連序75%的鹼基質量高於20

--max_barcode_errors barcode 容許的鹼基錯配個數，建議設置0爲不容許，即便修改成1，2一般也不容許錯配，不信你試試

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 # 全部序列長度分佈數據，可知長度範圍308-488,峯值爲408

seqs.fna # 質控並拆分後的數據，序列按樣品編號爲SampleID_0/1/2/3

split_library_log.txt # 日誌文件，有基本統計信息和每一個樣品的數據量；查看可知每一個樣品最大數據量爲110454，最小值爲10189

更多說明建議閱讀幫助http://qiime.org/scripts/split_libraries_fastq.html

 

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 --user

cutadapt切除雙端引物及長度控制，參數詳解：

-g 5’端引物

-a 3’端引物，須要將引物取反向互補

-e 引物匹配容許錯誤率，我調置0.15，通常引物20bp長容許3個錯配，爲了儘可能把引物切乾淨

-m 最小序列長度，根據狀況設置，本實驗擴增V5-V7區，長度主要位於350-400，故去除長度小於300bp的序列，無經驗可不填此參數

--discard-untrimmed 引物未切掉的序列直接丟棄，防止出現假OTU

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 3.6.1
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 73.83 s (58 us/read; 1.04 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.