ubiome类似数据dada2处理探索2

还是接着之前的事说，首先，在researchgte网站上发现了一个小“新大陆”，说可以把不能很好拼接的数据直接N连接处理。这里就先按软件默认的加几个N了，虽然拼接率有3%。。。然后，找到了直接加N的软件，不重复造轮子，自己写拼接脚本还是要费半天时间的，不如用工具，既好又准。这篇文章里的第一个软件死活找不到，于是放弃，还好有第二个选择，使用vsearch，依然是不能运行32程序的MacOS 10.15 Catalina，还是老样子，docker运行linux版本，至于为什么用32位，64位收费嘛。下面是我的操作过程：

1.“盲拼”序列

首先把usearch申请下载到工作目录，然后docker挂载到home，当然如果linux就直接省了这一步了，可以下载（安装）好直接使用。

docker run -it -v ~/Biodata/16S/process-ubiome-16S/primer-cutted:/home ubuntu
#开始拼接
./usearch -fastq_join sIL10-2-S-4_1.fastq -reverse sIL10-2-S-4_2.fastq -join_padgap NNNNNNNN   -join_padgapq IIIIIIII  -fastqout join.fq
tree #看下结果
.
├── join.fq
├── merged.fastq
├── sIL10-2-S-4_1.fastq
├── sIL10-2-S-4_1.fastq.gz
├── sIL10-2-S-4_2.fastq
├── sIL10-2-S-4_2.fastq.gz
└── usearch
ls -lh
total 277M
-rw-r--r-- 1 root root 115M Dec 30 02:47 join.fq
-rw-r--r-- 1 root root 2.7M Dec 30 02:38 merged.fastq
-rw-r--r-- 1 root root  63M Dec 30 02:26 sIL10-2-S-4_1.fastq
-rw-r--r-- 1 root root 9.0M Dec 30 02:24 sIL10-2-S-4_1.fastq.gz
-rw-r--r-- 1 root root  63M Dec 30 02:47 sIL10-2-S-4_2.fastq
-rw-r--r-- 1 root root 7.8M Dec 30 02:24 sIL10-2-S-4_2.fastq.gz
-rwxr-xr-x 1 root root 2.3M Nov 25 06:07 usearch

2.进入dada2流程

dada2流程，有嵌入到qiime2中的dada2插件，和单独的R包，实现效果应该是一样的，那就两种方法都试一下吧，反正，我对R语言不怎么熟悉，顺便熟悉下。

1）qiime2-dada2流程

为了简单，直接把序列重名为下机的标准文件名，这样序列导入就不需要制作文件名列表了。运行过程中出现了报错，估计不可行，看大部分序列长度是268，估计插件对于N的处理很严格，所以不能通过。

#激活conda
source /Volumes/MacOS/Miniconda3/bin/activate
#激活工作qiime2环境
conda activate qiime2-2019.10
#重命名文件
gzip join.fq #或者 pigz join.fa 如果装了多线程gzip pigz
mkdir data && mv join.fq.gz data/join_S67_L001_R2_001.fastq.gz
#导入数据
qiime tools import \
--type 'SampleData[SequencesWithQuality]' \
--input-path data \
--input-format CasavaOneEightSingleLanePerSampleDirFmt \
--output-path demux-single-end.qza
### quality control #visualization 
qiime demux summarize \
 --i-data demux-single-end.qza\
 --o-visualization microbiome.qzv
 ##filter 聚类
qiime dada2 denoise-single \
  --i-demultiplexed-seqs demux-single-end.qza \
  --p-trim-left 0 \
  --p-trunc-len 268 \
  --o-representative-sequences rep-seqs-dada2.qza \
  --o-table table-dada2.qza \
  --o-denoising-stats stats-dada2.qza \
  --p-n-threads 4
  报错Plugin error from dada2:

  No reads passed the filter. trunc_len (250) may be longer than read lengths, or other arguments (such as max_ee or trunc_q) may be preventing reads from passing the filter.

Debug info has been saved to /var/folders/05/zk0fr2h91h7_qsnr1kn0r72m0000gn/T/qiime2-q2cli-err-f8b4l_mr.log

2）R包dada2流程

Yes, although we strongly recommend using primers that allow your paired-end reads to overlap significantly. The mergePairs(..., justConcatenate=TRUE) option allows the paired reads to be joined without any overlap, but with 10Ns inserted in between the forward and reverse reads. The chimera removal and assignTaxonomy functions will handle such merged reads, although some other functions may fail (eg. addSpecies).

参考自官方文档FAQ，应该说比较权威了。这个就不需要其他软件了，这个包直接以10个N完成连接序列，然后修改下默认参数就好了。那么走一遍官方步骤，虽然说物种注释可能会失败，试试啦。很幸运，基本上完成了分析，看来这个方法还是可行的，不失为一个方法，虽然官方不推荐，毕竟，变异区中间缺失了大概10bp的样子。以下内容基本上来自官方文档，更改了部分，适合跑一个样本和以10个N拼接序列。

Ps。基本上写到最后了才发现宏基因组公众号刘永鑫大哥已经把文档在今年二月份翻译过了，如果直接看他的翻译应该节约不少时间，应该事前搜索下防止多废功夫呀，当然，看看英语练习下也不错。特别是我的序列是那种不能拼接的情况，和官方示例数据还有一些不同的。

准备工作–安装包

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("dada2", version = "3.10")
报错提示需要升级，于是运行下面代码升级后继续
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor") #清华镜像速度杠杠的
BiocManager::install(version = '3.10') #完成后重新运行上面代码
library(dada2); packageVersion("dada2")
#[1] ‘1.14.0’
path <- "/Users/zd200572/Biodata/16S/process-ubiome-16S/primer-cutted"

不希望一篇太长了，于是在这里分开下，后面继续，省得翻起来烦，就这样。