盘点扩增子序列拼接工具和方法
整理：谢鹏昊
修改：文涛
∙序列拼接
▪一.QIIME
▪二. Vsearch
▪三. Usearch
序列拼接
拿到扩增子测序结果之后的第一件事就是序列拼接（默认公司发放的数据都是没有问题的，所以不进行质量评估）。关于序列拼接，目前我们使用很多工具，三大扩增子分析工具（Qiime，mothur，usearch（vsearch））都有自己的拼接命令。目前mothur出现在文章中的
频率不高了，早先的软件跟新不够，相比于usearch不够快，相比于Qiime不够全面。所以，目前文章中做扩增子分析常用的也就是Qiime和usearch的拼接命令。Vsearch是一个怎样的存在呢？其实第一次听说，实在17年的nuture中，发布全球细菌微生物图谱。作为开源的扩增子分析工具，这个工具来源于usearch，功能和命令都类似。价值主要在可以免费使用，无数据量的限制。
所以下面就着三个工具 Qiime usearch Vsearch 的序列拼接命令进行一个总结，我们提供两种模式的序列拼接，基于单个fq文件，基于多个fq文件。（末尾有彩蛋）
一.QIIME
1.单个fq文件拼接
join_paired_ends.py -m SeqPrep –f $PWD/forward_reads.fastq -r $PWD/reverse_reads.fastq -o $PWD/SeqPrep_joined
join_paired_ends.py QIIME 软件上一个函数，执行拼接命令
-m 调用方法
1.SeqPrep 慢 更准确
2.fastq-join 快
–f 设置序列前端所在文件路径
-r 设置序列后端所在文件路径
-o 合并后输出的文件位置
2.多条序列同时拼接
# 对文件名有固定要求ABC_L001_R1_001 ABC_L001_R2_001
multiple_join_paired_ends.py -i input_files -o output_folder -p --include_input_dir_path
-i 输入文件目录 -o 输出文件目录
∙ERROR1
Invalid filename found for splitting on input for file /home/qiime/Desktop/Shared_Folder/NCBI28/16S/data_processing_gg135/a0_raw_data/a1_name_change/SRR9620715_1.fq, check input read1_indicator and read2_indicator parameters as well.
通过修改序列名字为标准格式解决
∙ERROR2Text file busy:'/home/qiime/Desktop/Shared_Folder/NCBI28/16S/data_processing_gg135/a0_raw_data/a1_name_change/try/S0_L001_R1_001/fastqjoin.un2
将文件夹try移动到桌面解决
3. for循环+单个拼接=多个fq文件拼接
for i in *_1.fq; do echo ${i%%_*}; join_paired_ends.py -m SeqPrep -f $PWD/ts 数组字符串转数组${i%%_*}_1.fq -r $PWD/${i%%_*}_2.fq -o $PWD/resistant_join1/${i%%_*}; mv $PWD/resistant_join1/${i%%_*}/seqprep_ $PWD/resistant_join2/${i%%_*}.; gzip -d $PWD/resistant_join2/${i%%_*}.; done
首先设置好工作目录
for i in; do ;done对目录中文件依次执行do 后面所列命令
for i in *_1.fq 将目录下所有名称为_1.fq 赋值给i
${i%%_*} 移除第一个_及其右边字符
eg: i= SRR9620715_1.fq --- SRR9620715
echo 在面板中显示结果
上面这条对多个fq文件拼接的命令拆分展示
#序列拼接
join_paired_ends.py -m SeqPrep -f $PWD/${i%%_*}_1.fq -r $PWD/${i%%_*}_2.fq -o $PWD/resistant_join1/${i%%_*}
#将合并后文件移动到另一文件夹，方便后续操作 mv a(文件所在) b(指定文件夹)
mv $PWD/resistant_join1/${i%%_*}/seqprep_ $PWD/resistant_join2/${i%%_*}.
# 对gz 文件解压 done 结束循环
gzip -d $PWD/resistant_join2/${i%%_*}.; done
二. Vsearch
1.单个拼接
# \为换行符,否则分行后无法运行
vsearch --fastq_mergepairs seq/WT1_1.fq \
--reverse seq/WT1_2.fq \
--fastqout d.fq \
--relabel WT1.
vsearch --fastq_mergepairs 执行拼接
seq/WT1_1.fq 序列前端所在文件位置
--reverse seq/WT1_2.fq 序列后端所在位置
--fastqout d.fq 输出目录
--relabel WT1. 改名
2. for循环+单个命令=多条拼接
这个循环需要我们提前准备好metadata文件，就是fq文件的名称列表。
for i in `tail -n+2 metadata.tsv | cut -f 1`;do
vsearch --fastq_mergepairs seq/${i}_1.fq --reverse seq/${i}_2.fq \
--fastqout temp/${i}.merged.fq --relabel ${i}.
done
tail -n+2去表头，cut -f 1取第一列，即获得样本列表
vsearch --fastq_mergepairs 拼接命令
seq/${i}_1.fq 序列前端
--reverse seq/${i}_2.fq 序列后端
--fastqout temp/${i}.merged.fq 输出路径及名字
--relabel ${i}.
三. Usearch
1.单个拼接
usearch -fastq_mergepairs SampleA_R1.fastq -reverse SampleA_R2.fastq -fastqout merged.fq
2.多个拼接
这个命令出了将序列拼接完成之后顺便讲全部拼接好的结果合并为一个gq文件，方便后续处理。
for Sample in *_R1.fq; do echo ${Sample%_*}; usearch -fastq_mergepairs ${Sample%_*}_R1.fq -fastqout $Sample.merged.fq -relabel $Sample.; cat $Sample.merged.fq >> d.fq;done
新一代的扩增子流程
基于新一代扩增子流程，Qiime2已经将序列拼接命令合并到的一起了，目的就是寻求最为简易封装的命令模式。我们可能也碰不到其中的序列拼接命令。其实这其中有一些细节需要注意，尤其是对于DADA2算法，这个算法并不是使用的传统的流程：序列拼接-质控-otu表格-下游的模式，而是使用：序列质控-DADA2错误学习-unique序列获得-序列拼接-去除嵌合的顺序做的，所以有很多朋友问询，是否可以将序列拼接后的结果直接做DADA2呢，这里不建议这么做。
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论