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

2021-03-03 宏基因組

點擊上方藍色「宏基因組」關注我們！專業乾貨每日推送！

寫在前面

之前發布的《擴增子圖表解讀》系列，相信關注過我的朋友大部分都看過了(連結直達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 密碼：y33d。

本課程代碼的運行，至少需要Linux平臺+安裝QIIME1.9.1，我之前發布過三種安裝QIIME的方法詳見文章目錄，總有一款適合你。

第二節. 提取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列，不用考慮方向。

# 提取barcodeextract_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%準確
—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,峰值為408seqs.fna # 質控並拆分後的數據，序列按樣品編號為SampleID_0/1/2/3split_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 2.7.12Command line parameters: -g AACMGGATTAGATACCCKG -a GGAAGGTGGGGATGACGT -e 0.15 -m 300 --discard-untrimmed temp/PE250_split/seqs.fna -o temp/PE250_P5.faTrimming 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,436Reads 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 bpTotal 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: 2Bases 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 sequenceslength    count    expect    max.err    error counts3    3    19959.9    0    34    4    4990.0    0    46    2    311.9    0    218    202626    0.0    2    20262619    56    0.0    2    56=== 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: 2Overview of removed sequenceslength    count    expect    max.err    error counts18    327    0.0    2    189 117 2119    1059175    0.0    2    1056863 2069 24320    13846    0.0    2    1817 10955 107421    744    0.0    2    5 10 729WARNING:    One or more of your adapter sequences may be incomplete.    Please see the detailed output above.

寫在後面

今天先到這裡，本文已經講了三個程序的使用，夠大家學習一會的了。要想了解這些程序的更多功能，一定要閱讀程序的幫助全文，才能有更深入的理解。

下節預告：擴增子分析解讀3去冗餘,聚類,生成OTU表

(宏基因組7月文章目錄，更多精彩等你讀)

Reference

http://qiime.org/scripts/index.html

http://qiime.org/scripts/extract_barcodes.html

http://qiime.org/scripts/split_libraries_fastq.html

https://pypi.python.org/pypi/cutadapt

想了解更多宏基因組、16S分析相關文章，

快關注「宏基因組」公眾號，乾貨第一時間推送。

系統學習生物信息，快關注「生信寶典」，

那裡有幾千志同道合的小夥伴一起學習。

相關焦點

  • 擴增子分析解讀4去嵌合體,非細菌序列,生成代表性序列和OTU表
    寫在前面之前發布的《擴增子圖表解讀》系列，相信很多朋友都看過了。這些內容的初衷是寫給剛進實驗室的學生讀，加速大家對同行文章的解讀能力。如果連同行的結果都看不懂，何談對數據的理解，對科學問題的解釋。希望剛入行的小夥伴多讀高水平文章，配合我的解讀，定能讓理解上升一個層次。
  • 16S擴增子測序數據分析介紹
    https://v.qq.com/x/page/t3015tp7d5u.htmlPart 1. 21擴增子分析背景介紹p1-23，23minhttps://v.qq.com/x/page/j3015gkf92g.htmlPart 2. 21擴增子分析背景介紹p24-39，32minhttps://v.qq.com/x/page/m3015en6o80
  • 微生物組-擴增子16S分析和可視化(線上/線下同時開課,2021.4)
    本課程一共3天，每天6節課，共18節課，全部課程均理論與實戰結合(只要課上講的內容，都是要帶你親自實現的分析)。從分析平臺搭建、Linux和R基礎、圖表解讀和繪圖實戰、擴增子分析標準流程、功能預測、差異統計分析以及各類高級分析(進化樹、網絡、環境因子、隨機森林、Adaboost和來源追溯等)，和CNS級圖片編輯和排版。
  • 微生物組-擴增子16S分析和可視化第10期(線上/線下同時開課,本年最後一期)
    本課程一共3天，每天6節課，共18節課，全部課程均理論與實戰結合(只要課上講的內容，都是要帶你親自實現的分析)。從分析平臺搭建、Linux和R基礎、圖表解讀和繪圖實戰、擴增子分析標準流程、功能預測、差異統計分析以及各類高級分析(進化樹、網絡、環境因子、隨機森林、Adaboost和來源追溯等)，和CNS級圖片編輯和排版。
  • 【選修一導學】 5.2 多聚酶鏈式反應擴增DNA片段
    2.提示：PCR引物是根據需要擴增的目標DNA的鹼基序列來設計的。(4)與體內DNA複製相比,PCR反應體系除引物外,還需加入哪種特別的物質?　　　　　　　　。 14.在遺傳病及刑偵破案中常需要對樣品DNA進行分析,PCR技術能快速擴增DNA片段。圖L5-2-5為某一次循環的產物,請據圖回答問題:
  • 擴增子圖表解讀2散點圖:組間整體差異分析(Beta多樣性)
    背景介紹(Introduction)宏基因組學宏基因組學目前的主要研究方法包括：16S/ITS/18S擴增子、宏基因組、宏轉錄組和代謝組，其中以擴增子研究最為廣泛。將來在大家可以很好理解相關文章圖表的基礎上，希望對分析、統計和繪圖相關技術有進一步學習的小夥伴請積極回復並留言吧。如果本系統文章閱讀過萬，想學分析的留言過百。我還將詳細講解擴增子分析、統計和繪圖各步驟的分析實例和原始碼，希望大家多多鼓勵和支持。聲明：文章的解讀僅代表個人理解和觀點，有不足處，請讀者積極留言批評指正，互相學習，共同進步。
  • 擴增子分析QIIME2. 1簡介和安裝
    擴增子分析QIIME2. 1簡介和安裝QIIME2版本 2017.6聲明：本文為QIIME2官方幫助文檔的中文版，
  • 相對豐度會歪曲實際豐度,聯合16S擴增子測序和總菌qPCR獲得的絕對豐度可靠嗎?
    發表時間：2020年發表單位：華盛頓大學弗雷德哈欽森癌症研究中心疫苗和傳染病部儘管16S擴增子測序可定量分析細菌類群的相對豐度，但樣品之間總細菌負荷的差異限制了其反映單個細菌物種絕對豐度的能力。qPCR可以定量單個物種的絕對豐度，但是針對存在於不同樣品中的每種細菌開發一套qPCR檢測方法是不切實際的。
  • 16S擴增子分析專題研討論會——背景介紹
    https://v.qq.com/x/page/t3015tp7d5u.htmlPart 1. 21擴增子分析背景介紹p1-23，23minhttps://v.qq.com/x/page/j3015gkf92g.htmlPart 2. 21擴增子分析背景介紹p24-39，32minhttps://v.qq.com/x/page/m3015en6o80
  • 「科學盒子」課程系列2|DNA的粗提取和PCR的擴增技術
    我們已經輕鬆地把DNA粗提取的問題解決了，下來就要研究名字很奇怪的「PCR」了。 PCR一般指聚合酶鏈式反應，所以可以簡單理解為PCR擴增技術就是把一小部分DNA無線擴增的技術。有了這項技術我們就可以把僅有的DNA擴增，例如把恐龍的DNA擴增，哇塞，我的寵物可不可以是恐龍呢？想想都很激動呢。
  • 圖解利用PCR技術擴增目的基因的過程
    在人教版高中生物教材選修一《專題5 DNA和蛋白質技術》和選修三《專題1 基因工程》中都介紹了這項技術，該知識點在2018年高考考綱中被列為選考部分II級要求，即能在複雜情境中綜合運用其進行分析、判斷、推理和評價。在新課講解時，該知識點的原理容易被學生理解，但涉及到該知識點的綜合運用型問題，一些同學暴露出了對細節把握不準的問題。
  • Realtime PCR數據分析注意事項以及數學原理
    1．理論上來說，應該使用目的基因和內參基因擴增至少5個一般7個濃度梯度稀釋的同一模板，如果兩者的擴增效率相等，那麼以目的基因濃度對ΔCT作圖，該直線的斜率理論值應該為0，實際值應該無限接近於0；2.第一，最重要的，溶解曲線，必須保證你的溶解曲線是單峰曲線，有任何一個小峰都不行，如果出現在主峰之前，考慮為引物二聚體，如果出現在主峰之後，考慮非特異性擴增，如果曲線完全沒有主峰，考慮引物汙染；第二，擴增曲線，首先必須保證你的syber green信號基本從0開始，如果不從0開始，考慮為基因組
  • 「實驗」PCR擴增目的片段
    它不僅可用於基因分離、克隆和核酸序列分析等基礎研究，還可用於疾病的診斷。PCR基本原理：類似於DNA的體內複製。在模板DNA、引物和4種脫氧核苷酸存在條件下依賴於DNA聚合酶的酶促合成反應。2、退火(Annealing)：當溫度突然降低時，反應體系中引物和其互補的DNA模板在局部形成雜交鏈。
  • 重大突破|日本科學家建立新式PCR擴增方法,能在常溫37°C快速擴增DNA
    DNA解離，使之成為單鏈，以便它與引物結合，為下輪反應做準備； 模板DNA與引物的退火(復性)：模板DNA經加熱變性成單鏈後，溫度降至55℃左右，引物與模板DNA單鏈的互補序列配對結合； 引物的延伸：DNA模板--引物結合物在Taq酶的作用下，以dNTP為反應原料，靶序列為模板，按鹼基配對與半保留複製原理，合成一條新的與模板DNA
  • DADA2中文教程v1.8
    開始之前本流程假設您的測序數據符合以下規範：樣品已被拆分好，即每個樣品一個fq/fastq文件（或者雙端成對fq文件）；已經去除非生物核酸序列，比如：引物（primers），接頭（adapters or barcodes），linker等；如果樣品是下機的雙端測序，其應具有雙端測序的相匹配的兩個fq文件。
  • 反義引物應該設置在引發cDNA合成的寡核苷酸的上遊
    ●cDNA 合成可用一種合成的反 義寡核苷酸引物來引發，該引物可選擇與特殊靶RNA或mRNA家族序列的某個區域雜交。特異cDNA片段能通過PCR擴增來獲得，而正義和反義寡核苷酸引物可根據cDNA具體序列來設計。
  • 微生物擴增子測序圖表解讀(實例數據)
    7.2 PCA分析主成分分析（Principal component analysis）PCA 是一種研究數據相似性或差異性的可視化方法，通過一系列的特徵值和特徵向量進行排序後，選擇主要的前幾位特徵值，採取降維的思想，PCA 可以找到距離矩陣中最主要的坐標，結果是數據矩陣的一個旋轉，它沒有改變樣品點之間的相互位置關係，只是改變了坐標系統。