ubiome类似数据dada2处理探索3

继续之前未完成的笔记,前面实现了使用qiime2-dada2插件初步探索,结果以无敌的报错失败告终,这里进入R包,更灵活地处理数据,下面是我的详细步骤。

1.上次未完成的准备工作

#设置R包的工作目录,不是R语言的,这只是个变量而已
path <- "/Users/zd200572/Biodata/16S/process-ubiome-16S/primer-cutted"
# 获得目录中的数据列表
fnFs <- sort(list.files(path, pattern="_1.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_2.fastq", full.names = TRUE))
# 获取样本名称
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
#看下样本,为了简单,我就放了一个样本
fnFs
[1] "/Users/zd200572/Biodata/16S/process-ubiome-16S/primer-cutted/sIL10-2-S-4_1.fastq"

2.数据质量质控

质量控制是必由之路,看下这个质量图,与qiime2的交互动态图不同,这个是静态的,

#看看质量情况,可以看出最后的质量值差到不忍直视了,难怪拼接率只有3%
plotQualityProfile(fnFs[1:2]) #正向质量
plotQualityProfile(fnRs[1:2]) #反向质量

3.过滤和修剪

这里还有个严重不合格样本的数据过滤,为了流程命名不变,我也走一遍。

# Place filtered files in filtered/ subdirectory
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
names(filtFs) <- sample.names
names(filtRs) <- sample.names

然后,我又通过fastqc看了下,序列基本上都是在120bp以后质量显著下降,于是正反向从120bp开始截取。

out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(120,120),
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=TRUE) # On Windows set multithread=FALSE
head(out)
                    reads.in reads.out
sIL10-2-S-4_1.fastq   198387    194949
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)
plotErrors(errF, nominalQ=TRUE)

这里说明一下,如果用之前我说的github上ubiome里面的数据的话,可能还要按照具体情况具体分析,应该是已经质控好的,直接加N即可,不过估计切到120再做也可以,就是损失几个碱基而已。

4.训练错误率模型?

因为我是一个样本,速度还是挺快的,几乎是秒完成。i3-3120M处理器,比较老啦,凑和还能用两年。然后做下来结果拟合度对比官方的感觉还挺好。

这是我的结果

官方的,区别不大吧

官方的,区别不大吧,所以应该是可以用的,直接继续了。

5.样本自参考

这应该是把上面建好的模型应用于我们的序列,这一步还是比较快的,也可能我只有一个样本的原因。作者提示如果windows系统要把multithread设置为FALSE。

dadaFs <- dada(filtFs, err=errF, multithread=TRUE)
#Sample 1 - 194949 reads in 11507 unique sequences.
dadaRs <- dada(filtRs, err=errR, multithread=TRUE)
#Sample 1 - 194949 reads in 11385 unique sequences.
dadaFs[[1]]
TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCGAGTTAAGTCAGCGGTAAAAGCCCGGGGCTCAACCCCGGCCCGCCGTTGAAACTGGCTGGCTTG                                                                                        35756 

6.合并双向reads

这里就是最最关键合并步骤了,justConcatenate=TRUE把序列直接以10个N拼接,亮眼的10个N。

mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE, justConcatenate=TRUE)
head(mergers[[1]]) #这里可以很明显地看到10个N
$TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCGAGTTAAGTCAGCGGTAAAAGCCCGGGGCTCAACCCCGGCCCGCCGTTGAAACTGGCTGGCTTG
[1] "TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCGAGTTAAGTCAGCGGTAAAAGCCCGGGGCTCAACCCCGGCCCGCCGTTGAAACTGGCTGGCTTGNNNNNNNNNNGCAGGCGGAATGCGCGGTGTAGCGGTGAAATGCATAGATATCGCGCAGAACCCCGATTGCGAAGGCAGCCTGCCGGCCCCACACTGACGCTGAGGCACGAAAGCGTGGGTATCGAACAGG"

7.构建序列表

这就是当年的OTU表,现在流行的ASV featuretable了,这里与之前OTU的区别应该是相似度,现在都是以100%相似度聚类了。

seqtab <- makeSequenceTable(mergers)
dim(seqtab)
[1]   1 500

8.去嵌合

去嵌合步骤中,作者提示如果大部分样本被除去,或者序列长度不对,要运行一下注释的那句。

table(nchar(getSequences(seqtab)))
#seqtab2 <- seqtab[,nchar(colnames(seqtab)) %in% 250:256]
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
head(seqtab.nochim)
#   TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCGAGTTAAGTCAGCGGTAAAAGCCCGGGGCTCAACCCCGGCCCGCCGTTGAAACTGGCTGGCTTGNNNNNNNNNNGCAGGCGGAATGCGCGGTGTAGCGGTGAAATGCATAGATATCGCGCAGAACCCCGATTGCGAAGGCAGCCTGCCGGCCCCACACTGACGCTGAGGCACGAAAGCGTGGGTATCGAACAGG
#[1,]                                                                                                                                                                                                                                                      35670
sum(seqtab.nochim)/sum(seqtab)
[1] 0.9975397

9.最后的检查

这里相当于汇总一下序列从原始到各个步骤,去除了多少条。

getN <- function(x) sum(getUniques(x))
track <- cbind(out, getN(dadaFs), getN(dadaRs), getN(mergers), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
head(track)

10.物种注释

#下载数据库,第一次运行
#download.file("https://zenodo.org/record/158955/files/gg_13_8_train_set_97.fa.gz?download=1",destfile = 'gg_13_8_train_set_97.fa.gz')
taxa <- assignTaxonomy(seqtab.nochim, "gg_13_8_train_set_97.fa.gz", multithread=TRUE)
#这个步骤不好用,可能是N太多导致不能更好分析到种,应该说对于这种盲目拼接的序列,到属是最好的选择,因为可能存在序列+10N之后拼接起来indel有多个。还可能是greengenes本身这个训练集不能到种,需要换另外的RDP或者Silva,以后再处理这个问题,而且后面我选择使用qiime2进行物种注释,发现这样结果更好些。
#taxas <- addSpecies(taxa, "/Users/zd200572/Biodata/16S/process-ubiome-16S/primer-cutted/rdp_species_assignment_16.fa.gz")
#Error in .Call2("ACtree2_build", tb, pp_exclude, base_codes, nodebuf_ptr,  : 
#  non base DNA letter found in Trusted Band for pattern 1
taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)[1,]
           Kingdom             Phylum              Class              Order 
     "k__Bacteria" "p__Bacteroidetes"   "c__Bacteroidia" "o__Bacteroidales" 
            Family              Genus            Species 
        "f__S24-7"              "g__"              "s__" 
#虽然这个只注释到科,还是可以接受的,毕竟v4只有1/8 16S全长序列。

12.导出序列

因为我R语言水平不咋的,就直接使用刘大哥的代码了:

#保存导出文件
setwd(path)
seqtable.taxa.plus <- cbind('#seq'=rownames(taxa), t(seqtab.nochim), taxa)
# ASV表格导出
write.table(t(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)
#————————————————
#版权声明:本文为CSDN博主「刘永鑫Adam」的原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接及本声明。
#原文链接:https://blog.csdn.net/woodcorpse/article/details/86773312

发表评论