利用dada2处理16S数据

生信喵

一、软件介绍
        DADA2(Divisive Amplicon Denoising Algorithm 2)是一个用于建模和修正多种测序平台（Illumina、Roche 454）测序扩增子错误的开源软件（R）包。在扩增子分析处理流程中，
        DADA2 算法能够准确地推断样本序列并寻找出单个核苷酸的差异(往往能够比其他方法识别出更多真实变体和输出更少的虚假序列）。
        DADA2 是一个通过构建错误率模型来推测扩增子序列是否来自模板的算法，以自身数据的错误模型为参数，不用依赖于其他参数分布模型。DADA2 最重要的优势是用了更多的数据，错误模型包含了质量信息，而其他的方法都在过滤低质量之后把序列的质量信息忽略。
        DADA2 的错误模型也包括了定量的丰度，而且该模型也计算了各种不同转置的概率。和比较 OTU 数据库的聚类算法不同，DADA2 采取的是降噪算法。聚类算法与降噪算法的差异与优劣。

1. 教程：http://benjjneb.github.io/dada2/tutorial.html
   
2. github：https://benjjneb.github.io/dada2/index.html
   
3. 案例：https://astrobiomike.github.io/amplicon/workflow_ex#analysis-in-r
   
4. 案例数据：https://mothur.s3.us-east-2.amazonaws.com/wiki/miseqsopdata.zip
   
5. 练习数据：
   
6. https://s3-eu-west-1.amazonaws.com/pfigshare-u-files/28773936/dada2_amplicon_ex_workflow.tar.gz

复制代码

dada2 要求数据具备下面三个条件。
        1、样品已被拆分好，即每个样品一个 fq/fastq 文件（或者双端成对 fq 文件）；
        2、已经去除非生物核酸序列，比如：引物（primers），接头（adapters or barcodes），linker等；
        3、如果样品是下机的双端测序，其应具有双端测序的相匹配的两个 fq 文件。


二、使用 dada2 分析 16S 数据
1、准备工作
        安装软件包

1. # install.packages('BiocManager',destdir = 'D:/R/Rstudiowork/downloaded_packages/')
   
2. #BiocManager::install("dada2",destdir = 'D:/R/Rstudiowork/downloaded_packages/')
   
3. #BiocManager::install("phyloseq",destdir = 'D:/R/Rstudiowork/downloaded_packages/')
   
4. #BiocManager::install("DECIPHER",destdir = 'D:/R/Rstudiowork/downloaded_packages/')
   
5. #BiocManager::install("vegan")
   
6. #BiocManager::install("msa",destdir = 'D:/R/Rstudiowork/downloaded_packages/')
   
7. # install.packages('shiny',destdir = 'D:/R/Rstudiowork/downloaded_packages/')
   
8. 
   
9. shiny::runGitHub("shiny-phyloseq","joey711")
   
10. install.packages("phangorn",destdir = 'D:/R/Rstudiowork/downloaded_packages/')
    

复制代码

       下载数据

1. # MiSeq SOP练习数据
   
2. #https://mothur.s3.us-east-2.amazonaws.com/wiki/miseqsopdata.zip
   
3. #DECIPHER注释库
   
4. #http://www2.decipher.codes/Classification/TrainingSets/SILVA_SSU_r138_2019.RData
   
5. 
   
6. #DADA2 SILVA库
   
7. #https://doi.org/10.5281/zenodo.4587954

复制代码

       加载软件包

1. library(dada2)
   
2. packageVersion("dada2")
   
3. library(phyloseq)
   
4. library(DECIPHER)
   
5. library(ggplot2)
   

复制代码

       设置工作目录

1. getwd()
   
2. path <- "./MiSeq_SOP"
   
3. list.files(path)

复制代码

2、读入数据

1. #定义reads1
   
2. fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
   
3. fnFs
   
4. #定义reads2
   
5. fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))
   
6. #提取样品名
   
7. sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
   
8. sample.names

复制代码

3、低质量过滤

1. #质控
   
2. plotQualityProfile(fnFs)
   
3. plotQualityProfile(fnRs)
   
4. 
   
5. #查看1-2样品
   
6. plotQualityProfile(fnFs[1:2])
   
7. 
   
8. #定义输出文件名
   
9. filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
   
10. filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
    
11. names(filtFs) <- sample.names
    
12. names(filtRs) <- sample.names
    
13. filtFs
    
14. filtRs
    
15. 
    
16. #过滤
    
17. out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
    
18.                      maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
    
19.                      compress=TRUE, multithread=FALSE)
    
20. # On Windows set multithread=FALSE
    
21. head(out)
    
22. 
    
23. #查看对象
    
24. class(out)
    
25. #查看维度
    
26. dim(out)
    
27. 
    
28. #过滤后质控
    
29. plotQualityProfile(filtFs)
    
30. plotQualityProfile(filtRs)
    
31. plotQualityProfile(filtFs[1:2])

复制代码

4、生成错误率模型

1. #生成reads1 错误率模型
   
2. errF <- learnErrors(filtFs, multithread=12)
   
3. #生成reads2 错误率模型
   
4. errR <- learnErrors(filtRs, multithread=12)
   
5. 
   
6. #可视化
   
7. plotErrors(errF, nominalQ=TRUE)

复制代码

5、降噪去重复

1. #核心步骤，降噪，聚类，生成ASV
   
2. dadaFs <- dada(filtFs, err=errF, multithread=12)
   
3. dadaRs <- dada(filtRs, err=errR, multithread=12)
   
4. 
   
5. dadaFs[[1]]
   
6. 
   
7. dadaFs$F3D0$clustering
   
8. str(dadaFs[[1]])

复制代码

6、合并 reads

1. mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)
   
2. #查看结果
   
3. head(mergers[[1]])
   
4. 
   
5. class(mergers)
   
6. length(mergers)
   
7. names(mergers)
   
8. 
   
9. class(mergers$F3D0)
   
10. names(mergers$F3D0)

复制代码

7、生成功能表

1. seqtab <- makeSequenceTable(mergers)
   
2. dim(seqtab)
   
3. #统计序列长度分布
   
4. table(nchar(getSequences(seqtab)))

复制代码

8、识别嵌合体

1. seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=12, verbose=TRUE)
   
2. dim(seqtab.nochim)

复制代码

9、统计结果

1. #定义函数
   
2. getN <- function(x) sum(getUniques(x))
   
3. track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
   
4. 
   
5. colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
   
6. rownames(track) <- sample.names
   
7. head(track)
   
8. 
   
9. #保存结果
   
10. write.table(track, "read-count-tracking.tsv", quote=FALSE, sep="\t", col.names=NA)

复制代码

憨憨1212

你好，我想我想请问一下。1、truncLen=c(240,160)是双端数据分别截取序列长度吗。这个和qiime2软件上对应的是--p-trunc-len-f和--p-trunc-len-r参数吗。2、还有就是truncLen=c(240,160)选择长度的依据是什么呢。加起来是400，所以是由于你的测序长度小于400bp的原因吗

生信喵

憨憨1212 发表于 2023-5-15 09:23
你好，我想我想请问一下。1、truncLen=c(240,160)是双端数据分别截取序列长度吗。这个和qiime2软件上对应的 ...

1对的，2选择的理由是前面质控的结果，R1超过240质量下降明显，R2超过160质量下降明显。