请选择 进入手机版 | 继续访问电脑版

生信人

找回密码
立即注册
搜索
热搜: 活动 交友 discuz
发新帖

0

收听

12

听众

246

主题
发表于 2023-5-11 17:16:27 | 查看: 1490| 回复: 3
一、软件介绍
        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. shiny::runGitHub("shiny-phyloseq","joey711")
  9. 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. #DADA2 SILVA库
  6. #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. #查看1-2样品
  5. plotQualityProfile(fnFs[1:2])

  6. #定义输出文件名
  7. filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
  8. filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
  9. names(filtFs) <- sample.names
  10. names(filtRs) <- sample.names
  11. filtFs
  12. filtRs

  13. #过滤
  14. out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
  15.                      maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
  16.                      compress=TRUE, multithread=FALSE)
  17. # On Windows set multithread=FALSE
  18. head(out)

  19. #查看对象
  20. class(out)
  21. #查看维度
  22. dim(out)

  23. #过滤后质控
  24. plotQualityProfile(filtFs)
  25. plotQualityProfile(filtRs)
  26. plotQualityProfile(filtFs[1:2])
复制代码

4、生成错误率模型
  1. #生成reads1 错误率模型
  2. errF <- learnErrors(filtFs, multithread=12)
  3. #生成reads2 错误率模型
  4. errR <- learnErrors(filtRs, multithread=12)

  5. #可视化
  6. plotErrors(errF, nominalQ=TRUE)
复制代码

5、降噪去重复
  1. #核心步骤,降噪,聚类,生成ASV
  2. dadaFs <- dada(filtFs, err=errF, multithread=12)
  3. dadaRs <- dada(filtRs, err=errR, multithread=12)

  4. dadaFs[[1]]

  5. dadaFs$F3D0$clustering
  6. str(dadaFs[[1]])
复制代码

6、合并 reads
  1. mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)
  2. #查看结果
  3. head(mergers[[1]])

  4. class(mergers)
  5. length(mergers)
  6. names(mergers)

  7. class(mergers$F3D0)
  8. 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. colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
  5. rownames(track) <- sample.names
  6. head(track)

  7. #保存结果
  8. write.table(track, "read-count-tracking.tsv", quote=FALSE, sep="\t", col.names=NA)
复制代码


收藏回复 显示全部楼层 道具 举报

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

回复 显示全部楼层 道具 举报

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

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

回复 显示全部楼层 道具 举报

您需要登录后才可以回帖 登录 | 立即注册

QQ|Archiver|手机版|小黑屋|生信人

GMT+8, 2024-3-28 23:29 , Processed in 0.132103 second(s), 20 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表