DADA2中文教程v1.8

EndNote相關資訊 | 2019-03-20 06:13

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

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

DADA2 Pipeline Tutorial

本文翻譯自DADA2英文使用指南:

我們將在一個小時內對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文件。

注意: 如果您的數據并不符合這些規范,那么需要您在開始本流程前解決這些問題,有關這類問題的建議,請參閱官方網站常見問題解答( ).

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

下面附上安裝代碼:

# 國內用戶清華鏡像站,加速國內用戶下載site=""# 檢查是否存在Biocondoctor安裝工具,沒有則安裝if (!requireNamespace("BiocManager", quietly = TRUE))  install.packages("BiocManager",repo=site)# 加載安裝工具,并安裝dada2library(BiocManager)BiocManager::install("dada2", version = "3.8")# 需要顯示版本信息library(dada2) packageVersion("dada2")備注:dada2的其他版本(ver. 1.2 1.4 1.8)均可同步實現該流程。

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

測試數據說明:這些腸道微生物測序樣品數據來自斷奶后的小鼠和一個是對照模擬群落(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中期望錯誤的最大值(參見“一種比平均質量更好的過濾策略”: );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"))# 過濾文件輸出,輸出和參數,統計結果保存于Outout <- 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.outF3D0_S188_L001_R1_001.fastq       7793      7113F3D1_S189_L001_R1_001.fastq       5869      5299F3D141_S207_L001_R1_001.fastq     5958      5463F3D142_S208_L001_R1_001.fastq     3183      2914F3D143_S209_L001_R1_001.fastq     3178      2941F3D144_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的流程參考: 目前鏈接無法打開,可能作者在更新)

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大數據工作流程( )。

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

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

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)# Inspect the merger data.frame from the first samplehead(mergers[[1]])mergers對象格式為R語言的數據框,數據框中包含序列及其豐度信息,未能成功合并的序列被刪除。

我們開始構建 amplicon sequence variant(ASV)表,;類似于我們傳統的OTU表一樣。

seqtab <- makeSequenceTable(mergers)dim(seqtab)# 查看序列長度分布# Inspect distribution of sequence lengthstable(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 nonchimF3D0    7793     7113      6976      6979   6540    6528F3D1    5869     5299      5227      5239   5028    5017F3D141  5958     5463      5331      5357   4986    4863F3D142  3183     2914      2799      2830   2595    2521F3D143  3178     2941      2822      2867   2552    2518F3D144  4827     4312      4151      4228   3627    3488序列物種注釋這里我們的16s序列注釋,采用下載Silva參考數據庫進行訓練和注釋,DADA2包為此目的提供了樸素貝葉斯分類器方法實現物種注釋。該assignTaxonomy函數將一組要分類的序列和具有已知分類的參考序列訓練集作為輸入,并輸出的注釋文件通過bootstrap檢驗。

首先從頁面 下載相應數據庫,有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 onlyrownames(taxa.print) <- NULLhead(taxa.print)不出所料,在這些糞便樣本中最豐富的是Bacteroidetes,但是僅僅有少數被注釋到種的水平,因為通常不可能從16S基因的一段進行確定的物種分配,并且可能因為參考數據庫中對小鼠腸道微生物群的覆蓋率極低。

另外:這里提到最近開發的IdTaxa物種注釋分類方法也可通過 DECIPHER Bioconductor包獲得。IDTAXA算法的論文表示其分類性能優于樸素貝葉斯分類器。這里我們包含一段代碼區域,允許你使用IdTaxa函數替代assignTaxonomy(并且它也更快?。?。經過訓練的分類器可從 獲得(可能需要翻墻)。下載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教程,并修改為下載目錄# 獲得(可能需要翻墻)load(paste0(path, "/SILVA_SSU_r132_March2018.RData")) # CHANGE TO THE PATH OF YOUR TRAINING SETids <- IdTaxa(dna, trainingSet, strand="top", processors=NULL, verbose=FALSE) # use all processorsranks <- 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 assignTaxonomytaxid <- 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 onlyrownames(taxa.print) <- NULLhead(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 Mockcat("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語言中有一整套的解決方案,盡請期待后續:++基于R語言phyloseq包的擴增子下游分析策略++。

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

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

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

系列教程:

專業技能:

一文讀懂:

必備技能:

文獻閱讀

擴增子分析:

在線工具:

科研經驗:

編程模板:

生物科普:

學習16S擴增子、宏基因組科研思路和分析實戰,關注“宏基因組”