搜索
查看: 512|回复: 0

ubiome数据探索(顺便测试图片md)

[复制链接]

8

主题

38

帖子

1009

积分

金牌会员

Rank: 6Rank: 6

积分
1009
QQ
发表于 2020-1-3 12:53:21 | 显示全部楼层 |阅读模式
  1. markdown

  2. 继续之前[未完成的笔记](https://mp.weixin.qq.com/s/Zw-4DasNB7ZBWmKZGQquow),前面实现了使用qiime2-dada2插件初步探索,结果以无敌的报错失败告终,这里进入R包,更灵活地处理数据,下面是我的详细步骤。

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

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

  16. #### 2.数据质量质控

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

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

  23. ![](https://github.com/zd200572/mp_kejijizhe/blob/master/blogs/Rplot01.png)

  24. ![](https://github.com/zd200572/mp_kejijizhe/blob/master/blogs/Rplot02.png)

  25. #### 3.过滤和修剪

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

  27. ```R
  28. # Place filtered files in filtered/ subdirectory
  29. filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
  30. filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
  31. names(filtFs) <- sample.names
  32. names(filtRs) <- sample.names
  33. ```

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

  35. ```R
  36. out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(120,120),
  37. maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
  38. compress=TRUE, multithread=TRUE) # On Windows set multithread=FALSE
  39. head(out)
  40. reads.in reads.out
  41. sIL10-2-S-4_1.fastq   198387    194949
  42. errF <- learnErrors(filtFs, multithread=TRUE)
  43. errR <- learnErrors(filtRs, multithread=TRUE)
  44. plotErrors(errF, nominalQ=TRUE)
  45. ```

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

  47. #### 4.训练错误率模型?

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

  49. ![这是我的结果](https://github.com/zd200572/mp_kejijizhe/raw/master/blogs/Rplot.png)

  50. ![官方的,区别不大吧](http://benjjneb.github.io/dada2/tutorial_files/figure-html/plot-errors-1.png)

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

  52. #### 5.样本自参考

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

  54. ```R
  55. dadaFs <- dada(filtFs, err=errF, multithread=TRUE)
  56. #Sample 1 - 194949 reads in 11507 unique sequences.
  57. dadaRs <- dada(filtRs, err=errR, multithread=TRUE)
  58. #Sample 1 - 194949 reads in 11385 unique sequences.
  59. dadaFs[[1]]
  60. TACGGAGGAGTTG...GAAACTGGCTGGCTTG         35756
  61. ```



  62. #### 6.合并双向reads

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

  64. ```R
  65. mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE, justConcatenate=TRUE)
  66. head(mergers[[1]]) #这里可以很明显地看到10个N
  67. $TACGGAGGATGCGAGCGTTATCCGG...CTCAACCCCGGCCCGCCGTTGAAACTGGCTGGCTTG
  68. [1] "TACGGAGGATGCGAGCGTTATCC...TGGCTGGCTTGNNNNNNNNNNGCAGGCGGA...GGTATCGAACAGG"
  69. ```

  70. #### 7.构建序列表

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

  72. ```R
  73. seqtab <- makeSequenceTable(mergers)
  74. dim(seqtab)
  75. [1]   1 500
  76. ```

  77. #### 8.去嵌合

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

  79. ```R
  80. table(nchar(getSequences(seqtab)))
  81. #seqtab2 <- seqtab[,nchar(colnames(seqtab)) %in% 250:256]
  82. seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
  83. head(seqtab.nochim)
  84. #   TACGGAGGATGCG...GGCTTGNNNNNNNNNNGCAGGCGG...TCGAACAGG
  85. #[1,]               35670
  86. sum(seqtab.nochim)/sum(seqtab)
  87. [1] 0.9975397
  88. ```



  89. #### 9.最后的检查

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

  91. ```R
  92. getN <- function(x) sum(getUniques(x))
  93. track <- cbind(out, getN(dadaFs), getN(dadaRs), getN(mergers), rowSums(seqtab.nochim))
  94. colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
  95. rownames(track) <- sample.names
  96. head(track)
  97. ```



  98. #### 10.物种注释

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

  116. #### 12.导出序列

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

  118. ```R
  119. #保存导出文件
  120. setwd(path)
  121. seqtable.taxa.plus <- cbind('#seq'=rownames(taxa), t(seqtab.nochim), taxa)
  122. # ASV表格导出
  123. write.table(t(seqtab.nochim), "dada2_counts.txt", sep="\t", quote=F, row.names = T)
  124. # 带注释文件的ASV表格导出
  125. write.table(seqtable.taxa.plus , "dada2_counts.taxon.species.txt", sep="\t", quote=F, row.names = F)
  126. # track文件保存
  127. write.table(track , "dada2_track.txt", sep="\t", quote=F, row.names = F)
  128. ```
复制代码



欢迎交流!
https://github.com/zd200572/
https://jiawen.zd200572.com
回复

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-7-14 01:02 , Processed in 0.033186 second(s), 31 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.