DADA2中文教程v1.8

DADA2 R包中文使用指南

譯者：文濤 南京農業的大學

責編：劉永鑫 中科院遺傳發育所

DADA2 Pipeline Tutorial

本文翻譯自DADA2英文使用指南：http://benjjneb.github.io/dada2/tutorial.html

寫在前面

我們將在一個小時內對DADA2（ver 1.8）包進行一個初步系統的學習。本教程使用的是已拆分好的雙端illumina測序下機數據，barcode或者接頭序列已經被去除。通過本流程我們將得到ASV 表格（amplicon sequence variant table），類似於我們熟悉的OTU table 。如同之前介紹過的DADA2特性，ASV table通過將每條序列及其數量信息記錄下來，得到具有比傳統聚類得到的OTU table更高解析度的微生物分類信息。通過注釋序列得到物種信息，下遊我們將使用目前十分流行的R包phyloseq做下遊微生物群落相關分析，盡請期待。

開始之前

本流程假設您的測序數據符合以下規範：

樣品已被拆分好，即每個樣品一個fq/fastq文件（或者雙端成對fq文件）；

已經去除非生物核酸序列，比如：引物（primers），接頭（adapters or barcodes），linker等；

如果樣品是下機的雙端測序，其應具有雙端測序的相匹配的兩個fq文件。

注意: 如果您的數據並不符合這些規範，那麼需要您在開始本流程前解決這些問題，有關這類問題的建議，請參閱官方網站常見問題解答（ http://benjjneb.github.io/dada2/faq.html ）.

準備

首先我們需要載入dada2包，如果您在此之前並沒有安裝成功，參見dada2安裝指南（http://benjjneb.github.io/dada2/dada-installation.html） 本來我準備跳過這個步驟，但是年底我一直在回家的途中，拿的是一臺新買的超極本給大家寫教程，機器上並沒有配置好dada2包，因此，這裡和沒有準備好這一環境的朋友們一起安裝。

下面附上安裝代碼:

# 國內用戶清華鏡像站，加速國內用戶下載
site="https://mirrors.tuna.tsinghua.edu.cn/CRAN"
# 檢查是否存在Biocondoctor安裝工具，沒有則安裝
if (!requireNamespace("BiocManager", quietly = TRUE))
 install.packages("BiocManager",repo=site)
# 加載安裝工具，並安裝dada2
library(BiocManager)
BiocManager::install("dada2", version = "3.8")

# 需要顯示版本信息
library(dada2)
packageVersion("dada2")

備註:dada2的其他版本（ver. 1.2 1.4 1.8）均可同步實現該流程。

數據獲取

官方地址：https://www.mothur.org/w/images/d/d6/MiSeqSOPData.zip

本翻譯教程全套流程文件亦包含所下載的原始數據

測試數據說明：這些腸道微生物測序樣品數據來自斷奶後的小鼠和一個是對照模擬群落（Illumina MiSeq，16S rRNA基因的2x250，V4區域），

定義原始文件路徑

定義工作路徑,顯示文件夾內容，這裡我將作者的linux路徑更換為win下的路徑；

# 設置數據所在目錄，請設置為你下載解壓的目錄
# CHANGE ME to the directory containing the fastq files after unzipping.
path <- "C:/Users/woodc/Desktop/home/markdown/blog/wentao/MiSeq_SOP"
list.files(path)

通過字符串操作函數提取單個樣品測序文件信息

# list.files返回指定目錄中的文件名
# pattern指定返回指定類型的文件
# full.names = TRUE返回帶有路徑的完整文件名
# 返回測序正向文件完整文件名
fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
# 返回測序反向文件完整文件名
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))

# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
# basename(fnFs)提取文件名(不要目錄)
# strsplit按`_`進行分割字符，返回列表
# sapply批量操作，這裡批量提取列表中第一個元素，即樣本名
# 提取文件名中`_`分隔的第一個單詞作為樣品名
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

序列文件質量檢測

繪製前4個樣本的質量示意圖

plotQualityProfile(fnFs[1:4])

這是一張熱圖，基於鹼基的質量信息的頻率分布所得，每個位置的中位數質量得分由綠線表示，質量得分的四分位數由橙色線表示。紅線顯示擴展至該位置序列長度的比例（這對於其他測序技術更有用，如454測序。因為Illumina讀數通常都是相同的長度，因此是扁平的紅線）。正向序列質量很好。我們通常建議修剪最後幾個核苷酸序列，最後幾個鹼基錯誤率往往比較高。我們將截取/保留前240位的正向序列（剪去最後10個核苷酸）。

反向序列的質量

現在我們可視化反向序列的質量概況

plotQualityProfile(fnRs[1:4])

反向序列質量明顯較差，特別是在最後，這在Illumina測序中很常見。但這對於我們的DADA2 流程影響不大。雖然DADA2將質量信息整合到其誤差模型中，但是這使得算法對低質量序列具有魯棒性，所以我們從平均質量值突降的位置進行剪切會提高算法對稀有序列的敏感性。因此，我們對反向序列文件裁剪去至160bp。

注意：當我們在處理自己的數據時，不僅僅需要根據質量突降的位置進行裁剪位置的選擇，首先我們必須確定測序的雙端序列必須擁有足夠的重疊(overlap)，至少保證裁剪之後在 20 bp 以上，以保證拼接效率。

序列過濾和裁剪

首先對序列進行質量過濾剪切。我們通過filterAndTrim命令對正向序列和反向序列分別裁剪至240和160 bp（truncLen=c(240,160)），最大N的數量設置為0（maxN=0），maxEE參數表示允許在條reads中期望錯誤的最大值（參見「一種比平均質量更好的過濾策略」：https://academic.oup.com/bioinformatics/article/31/21/3476/194979 ）；truncQ為過濾最最質量的閾值，默認為2，將去除低於此質量的序列；rm.phix默認為TRUE，將去除比對上參考PhiX基因組的序列，在擴增子中常用；compress，默認輸出壓縮格式的結果；multithread為默認單線程，可以改為T或整數提高速度，但確保你電腦有足夠的內存和計算資源。

# Place filtered files in filtered/ subdirectory
# 將過濾後的文件存於filtered子目錄，設置輸出文件名
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
# 過濾文件輸出，輸出和參數，統計結果保存於Out
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
                    maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                    compress=TRUE, multithread=T)
head(out)

顯示過濾前後的結果統計

                             reads.in reads.out
F3D0_S188_L001_R1_001.fastq       7793      7113
F3D1_S189_L001_R1_001.fastq       5869      5299
F3D141_S207_L001_R1_001.fastq     5958      5463
F3D142_S208_L001_R1_001.fastq     3183      2914
F3D143_S209_L001_R1_001.fastq     3178      2941
F3D144_S210_L001_R1_001.fastq     4827      4312

注意事項1：以上設置的過濾參數並不是不變的，如果我們需要縮短過濾時間，maxEE參數設置更小一些，想要保留多數的reads，那就需要對maxEE參數設置的更大些，尤其是反向測序數據（例如maxEE=c(2,5)）。另外，在設置trunclen參數時一定要注意在裁剪之後的序列能保留足夠的長度（最低20bp）來用於拼接。

注意事項2：對於ITS測序結果，序列長度變化較大，可以考慮不進行裁剪，但是要確保引物在這之前已經被去除乾淨。

錯誤率

下面我們來了解錯誤率

DADA2算法使用機器學習構建參數誤差模型，認為每個擴增子測序樣品都具有不同的誤差比率。該learnErrors方法通過交替估計錯誤率和對參考樣本序列學習錯誤模型，直到學習模型同真實錯誤率收斂於一致。與許多機器學習方法一樣，算法有一個初始假設：數據中的最大可能錯誤率就是只有最豐富的序列是正確的，其餘都是錯誤。（在我這臺小筆記本上運行用了6分鐘）

這是DADA2中運行最耗費計算資源的一步

errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)
plotErrors(errF, nominalQ=TRUE)
plotErrors(errR, nominalQ=TRUE)

可視化結果表示了每種可能的錯誤（A→C，A→G，…），圖中點表示觀察得到的錯誤率。黑線表示通過算法學習評估得到的錯誤率。紅色的曲線表示由Q-score的定義下預期的錯誤率。這裡估計的錯誤率（黑線）同觀察到的錯誤率（點）擬合程度很好，並且錯誤率隨著預期的質量而下降。這和我們的認知相符。

注意事項：DADA2核心算法亦是參數學習，計算量非常可觀，目前我們測得的數據大小為本例子文件的十倍甚至一百倍，面對如此巨大的數據和需要消耗的計算資源，這一模型的展示便不適合我們實際較大的數據量，可以通過增加nbase參數調整擬合程度以減少計算量。

去除重複序列

與usearch去冗餘步驟類似，僅僅保留重複序列中的一條序列，大量節省計算資源。對於在此使用小筆記本為大家翻譯教程的我十分受用。

在這裡值得注意的是DADA2保留了去重序列的質量信息，這些質量信息取自重複序列的均值。這一信息文件將作為參考錯誤模型用於後續序列處理，以提高了DADA2算法準確性。

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)

注意事項：如果您的數據量很大，可能需要使用別的策略更為妥當（big data使用DADA的流程參考：http://benjjneb.github.io/dada2/bigdata.html 目前連結無法打開，可能作者在更新）

基於錯誤模型進一步質控

dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
dadaFs[[1]]

DADA2算法從第一個樣本中的1979個獨特序列推斷出128個真實序列。dada-class返回對象還有許多（參見help(「dada-class」)參考資料），包括關於每個去噪序列質量的多個評價指標，本教程不做解釋。

注意事項1：在教程中，所有樣本都同時加載到內存中。如果您正在處理接近或超過可用內存（RAM）的數據集（如今我們的電腦內存普遍4或8G，據我測試，8G內存可處理18個樣品，每個樣品大小為20M，也就是兩個實驗組），樣品數量較多時，我們需要逐個處理樣本：請參閱 DADA2大數據工作流程（ http://benjjneb.github.io/dada2/bigdata.html ）。

注意事項2：DADA2同時支持454和IonTorrent數據，但是建議對其中一些參數進行修改，可以通過R語言中?setDadaOpt來調取幫助文件，探索這些參數的修改方法。（如有需求，歡迎留言交流）

序列拼接

我們現在將正向和反向序列文件合併在一起以獲得完整的序列。通過將去噪的正向序列與相應的去噪反向序列的反向互補序列比對，然後構建合併的「overlap」進行合併。默認情況下，僅當正向和反向序列重疊至少12個鹼基並且在重疊區域中彼此相同時才輸出合併序列。

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
# Inspect the merger data.frame from the first sample
head(mergers[[1]])

mergers對象格式為R語言的數據框，數據框中包含序列及其豐度信息，未能成功合併的序列被刪除。

生成ASV表

我們開始構建 amplicon sequence variant（ASV）表，；類似於我們傳統的OTU表一樣。

seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# 查看序列長度分布
# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

去除嵌合體

Dada核心質控方法糾正了一些錯誤，但嵌合體仍然存在。幸運的是，去噪後序列的準確性使得識別嵌合體比處理模糊OTU時更簡單。

seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochim)
sum(seqtab.nochim)/sum(seqtab)

嵌合體的數量歷來被大家討論過很多次，因為不同實驗，不同樣品等等很多的因素都於嵌合體數量有關，這裡我們得到的嵌合體數量大概為21%，但是這些嵌合體佔據的序列數量大部分都是很少的，在這裡這些嵌合體僅僅佔我們全部序列數量的3.6%。

注意事項: 在進行嵌合體去除操作後大部分序列應該被保留下來。也有可能大部分序列作為嵌合體被去除，但這樣情況並不多見。如果這種情況發生，我們只有返回上遊進行重新處理序列，大部分這種情況都是由於引物序列未能去除乾淨導致的。

統計上述分析步驟

以上步驟序列可算是進行了重重篩選，這些序列保留及其去除的數量信息需要做一個記錄才好。

# 求和函數
getN <- function(x) sum(getUniques(x))
# 合併各樣本分步數據量
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
# If processing a single sample, remove the sapply calls: e.g. replace sapply(dadaFs, getN) with getN(dadaFs)
# 修改列名
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
# 行名修改為樣本名
rownames(track) <- sample.names
# 統計結果預覽
head(track)

預覽前6個樣本各步過濾和處理前後數據量統計

      input filtered denoisedF denoisedR merged nonchim
F3D0    7793     7113      6976      6979   6540    6528
F3D1    5869     5299      5227      5239   5028    5017
F3D141  5958     5463      5331      5357   4986    4863
F3D142  3183     2914      2799      2830   2595    2521
F3D143  3178     2941      2822      2867   2552    2518
F3D144  4827     4312      4151      4228   3627    3488

序列物種注釋

這裡我們的16s序列注釋，採用下載Silva參考資料庫進行訓練和注釋，DADA2包為此目的提供了樸素貝葉斯分類器方法實現物種注釋。該assignTaxonomy函數將一組要分類的序列和具有已知分類的參考序列訓練集作為輸入，並輸出的注釋文件通過bootstrap檢驗。

首先從頁面 http://benjjneb.github.io/dada2/training.html 下載相應資料庫，有16S可選 Silva 132/128/123，RDP trainset 16/14, Greengene 13.8；真菌ITS必選UNITE。這裡我們使用Silva 132，需要下載序列訓練集和物種注釋兩個文件。

保存文件至工作目錄。其它位置請修改下面代碼中path變量

taxa <- assignTaxonomy(seqtab.nochim, paste0(path, "/silva_nr_v132_train_set.fa.gz"), multithread=TRUE)
taxa <- addSpecies(taxa, paste0(path, "/silva_species_assignment_v132.fa.gz"))
taxa.print <- taxa
# 另存物種注釋變量，去除序列名，只顯示物種信息
# Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)

不出所料，在這些糞便樣本中最豐富的是Bacteroidetes，但是僅僅有少數被注釋到種的水平，因為通常不可能從16S基因的一段進行確定的物種分配，並且可能因為參考資料庫中對小鼠腸道微生物群的覆蓋率極低。

另外：這裡提到最近開發的IdTaxa物種注釋分類方法也可通過 DECIPHER Bioconductor包獲得。IDTAXA算法的論文表示其分類性能優於樸素貝葉斯分類器。這裡我們包含一段代碼區域，允許你使用IdTaxa函數替代assignTaxonomy（並且它也更快！）。經過訓練的分類器可從 http://DECIPHER.codes/Downloads.html 獲得(可能需要翻牆)。下載SILVA SSU r132文件以繼續。（本教程未進行相關操作，各位有需求親留言交流）

# 安裝DECIPHER進行物種注釋
BiocManager::install("DECIPHER", version = "3.8")
library(DECIPHER); packageVersion("DECIPHER")

# 轉換ASV表為DNAString格式
dna <- DNAStringSet(getSequences(seqtab.nochim)) # Create a DNAStringSet from the ASVs
# 相關數據下載詳見DECIPHER教程，並修改為下載目錄
# http://DECIPHER.codes/Downloads.html 獲得(可能需要翻牆)
load(paste0(path, "/SILVA_SSU_r132_March2018.RData")) # CHANGE TO THE PATH OF YOUR TRAINING SET
ids <- IdTaxa(dna, trainingSet, strand="top", processors=NULL, verbose=FALSE) # use all processors
ranks <- c("domain", "phylum", "class", "order", "family", "genus", "species") # ranks of interest
# Convert the output object of class "Taxa" to a matrix analogous to the output from assignTaxonomy
taxid <- t(sapply(ids, function(x) {
 m <- match(ranks, x$rank)
 taxa <- x$taxon[m]
 taxa[startsWith(taxa, "unclassified_")] <- NA
 taxa
}))
colnames(taxid) <- ranks; rownames(taxid) <- getSequences(seqtab.nochim)
taxa.print <- taxa <- taxid# Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)
# 註：此處我沒有運行成功，讀入RData後R環境崩潰了，個人感覺R語言處理大文件非常不擅長

注意事項：如果您的數據沒有被適當注釋，例如您的細菌16S序列被分配為大量Eukaryota NA NA NA NA NA， 可能核苷酸序列方向與參考資料庫的方向相反。告訴dada2嘗試反向互補方向進行匹配，assignTaxonomy(…,tryRC=TRUE)看看這是否可以修復注釋信息。如果使用DECIPHER進行分類，請嘗試IdTaxa (…, strand=」both」)。

相關文件保存

setwd(path)
seqtable.taxa.plus <- cbind('#seq'=rownames(taxa), t(seqtab.nochim), taxa)
# ASV表格導出
write.table(seqtab.nochim, "dada2_counts.txt", sep="\t", quote=F, row.names = T)
# 帶注釋文件的ASV表格導出
write.table(seqtable.taxa.plus , "dada2_counts.taxon.species.txt", sep="\t", quote=F, row.names = F)
# track文件保存
write.table(track , "dada2_track.txt", sep="\t", quote=F, row.names = F)

DADA2準確性評估

本教程的最後我們來評估DADA2的準確性

在群落中我們包含了一個模擬的微生物群落，這是對目前已知的20株菌的混合測序樣品。並且已經提供了這些菌的真實注釋信息，我們使用DADA2進行推斷和真正的注釋信息進行比較，從而評估其注釋準確性。

這一人工群落包含20株菌DADA2確定了20個 ASV，所有這些ASV都與預期的群落成員的參考基因組完全匹配。此樣本的DADA2分析之後的殘留錯誤率為0％。

unqs.mock <- seqtab.nochim["Mock",]
unqs.mock <- sort(unqs.mock[unqs.mock>0], decreasing=TRUE) # Drop ASVs absent in the Mock
cat("DADA2 inferred", length(unqs.mock), "sample sequences present in the Mock community.\n")
mock.ref <- getSequences(file.path(path, "HMP_MOCK.v35.fasta"))
match.ref <- sum(sapply(names(unqs.mock), function(x) any(grepl(x, mock.ref))))
cat("Of those,", sum(match.ref), "were exact matches to the expected reference sequences.\n")

Dada2 piplines到此結束

我們得到了ASV表格，類似於傳統OTU 表格，注釋信息，傳統的代表序列文件在這裡的表現形式在ASV 表中的行名。我們除了沒有構建進化樹，基於擴增子的前處理就到此結束了。

基於R語言DADA2的後續分析

這裡我們在R語言中有一整套的解決方案，盡請期待後續：
++基於R語言phyloseq包的擴增子下遊分析策略++。

譯者簡介

個人簡歷：文濤，2016年就讀於南京農業的大學。在資源院沈其榮教授課題組，研究方向為根際微生物生態，具體為植物介導下根際小分子代謝組同土壤微生物群落在防控土傳病害方面的相互作用。2018年9月轉博，期間個人主要完成了一篇發表於Microbiome的文章內容（三作，一二做作為小導師及其合作關係）。題目為：Root exudates drive the soil-borne legacy of aboveground pathogen infection（https://doi.org/10.1186/s40168-018-0537-x） ；本人一作已投文章一篇，在寫文章一篇，兩篇文章主要數據均為代謝組和擴增子測序獲得。

技能：掌握擴增子測序傳統數據處理及其下遊分析；會使用Qiime2，usearch， R等工具進行新一代測序數據處理及其部分下遊分析。。想通過coding的學習或者大家的幫助將分析過程流程化。會基於R語言的非靶向代謝組學下遊分析技術。

現狀：目前在跟進碩士內容，並加入根際功能方面的內容（宏基因組），其中水稻相關的根際宏基因組樣品與12月送樣，希望熟練掌握宏基因組數據處理流程。同時希望完善代謝組學分析技術和基於擴增子和代謝組的大數據整合技術。

猜你喜歡

10000+：菌群分析 寶寶與貓狗 梅毒狂想曲 提DNA發Nature Cell專刊 腸道指揮大腦

系列教程：微生物組入門 Biostar 微生物組  宏基因組

專業技能：學術圖表 高分文章 生信寶典 不可或缺的人

一文讀懂：宏基因組 寄生蟲益處 進化樹

必備技能：提問 搜索  Endnote

文獻閱讀 熱心腸 SemanticScholar Geenmedical

擴增子分析：圖表解讀 分析流程 統計繪圖

16S功能預測   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在線工具：16S預測培養基 生信繪圖

科研經驗：雲筆記  雲協作 公眾號

編程模板: Shell  R Perl

生物科普:  腸道細菌 人體上的生命 生命大躍進  細胞暗戰 人體奧秘  

寫在後面

為鼓勵讀者交流、快速解決科研困難，我們建立了「宏基因組」專業討論群，目前己有國內外5000+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註「姓名-單位-研究方向-職稱/年級」。PI請明示身份，另有海內外微生物相關PI群供大佬合作交流。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。

學習16S擴增子、宏基因組科研思路和分析實戰，關注「宏基因組」

點擊閱讀原文，跳轉最新文章目錄閱讀