1.概览
# 目录
mkdir -p temp # 临时文件
mkdir -p result # 最终结果
# 文件
dp_16s_v18.fa #16S参考数据库第18版（可在http://www.drive5.com/sintax/rdp_16s_v18.fa.gz获得）
seq/*.fq.gz #压缩的原始测序数据
doc/design.txt #实验设计文件
2.PE reads拼接（fastq_mergepairs）
# 测序数据解压
gunzip seq/* 或 unpigz seq/*
# 依照实验设计批处理并合并
for i in `tail -n+2 doc/design.txt | cut -f 1`;
do
vsearch --fastq_mergepairs seq/${i}_1.fq --reverse seq/${i}_2.fq \
--fastqout temp/${i}.merged.fq \
--relabel ${i}. #序列重命名
done
# 合并所有样品至同一文件
cat temp/*.merged.fq > temp/all.fq
3.质量控制（fastx_filter）
# 去除引物（请按实际修改，如Cut barcode 10bp + V5 19bp in left and V7 18bp in right）
vsearch --fastx_filter temp/all.fq \
--fastq_stripleft 29 --fastq_stripright 18 \
--fastq_maxee_rate 0.01 \ #expected error #大于1.0则去除
--fastaout temp/filtered.fa
- 去冗余与聚类生成OTUs或降噪生成ASVs（
derep_fulllength，cluster_fast，uchime_ref，usearch_global）
#序列去冗余，推荐使用vsearch，并添加miniuniqusize为8，去除低丰度，增加计算速度
vsearch --derep_fulllength temp/filtered.fa \ #序列去冗余
--sizeout \ #输出结果带丰度信息
--minuniquesize 8 \ #过滤低丰度reads，常用8，RPM1
--output temp/uniques.fa #输出结果
#聚类方式生成OTUs或降噪生成ASVs
#按丰度高到低聚类选择cluster_fast，非聚类的精度序列变异选择cluster_unoise算法
vsearch --cluster_fast temp/uniques.fa \ #按序列长度排序聚类
--id 0.97 \ #相似性阈值
--centroids temp/otus.fa \ #输出中心序列作为代表序列
--relabel OTU_ #OTU序列重命名
vsearch cluster_unoise temp/uniques.fa \#采用unoise3算法去噪
--centroids temp/asvs.fa \ #输出中心序列作为代表序列
--relabel ASV_ #OTU序列重命名
#去除嵌合体
#细菌可用Usearch作者整理的RDP Gold数据库去除嵌合体（http://drive5.com/uchime/rdp_gold.fa）
vsearch --uchime_ref temp/otus.fa \ #基于参考数据库去嵌合
--db db/rdp_gold.fa \ #使用–uchime_ref 时指定数据库fasta文件
--nonchimeras result/otus.fa #输出无嵌合体结果文件
vsearch --uchime3_denovo temp/asvs.fa \#序列自身比对去嵌合
--abskew 16 #亲本丰度比是嵌合体16倍以上
--nonchimeras result/asvs.fa #输出无嵌合体结果文件
#物种注释，上面为OTU，下面为ASV（不建议，建议使用后面的方法）
vsearch --sintax otus.fa \
--db rdp_16s_v18.fa \ #前面下载的16S参考数据库
--tabbedout otutab_tax_anno.txt
vsearch --sintax asvs.fa \
--db rdp_16s_v18.fa \ #前面下载的16S参考数据库
--tabbedout asvtab_tax_anno.txt
5.统计OTUs或ASVs的计数值
#前面几步相当于建立一个与此实验或样本关联的专有OTUs数据库
vsearch --usearch_global temp/filtered.fa \ #全局比对之检索的序列
--db result/otus.fa \ #全局比对之数据库
--id 0.97 \ #相似性阈值：当查询序列与目标序列之间的相似度达到多少时，才算比对上
--query_cov 0.97 \ #覆盖度：满足相似度的情况下，同时要求查询序列的覆盖度达到多少；
--strand both \ #默认只检测正链，此处改为双链
--otutabout result/otutab.txt #经典表格格式 OTU表
推荐的物种注释方法（在统计完OTUs数目之后），物种注释使用参数与统计OTUs数目时一致（参考自https://zhuanlan.zhihu.com/p/242638228）
vsearch --usearch_global ${outdir}/${sample}_otu.fa \
-db $DB/amplicon.fa \ #16S参考数据库
--id 0.97 \
--query_cov 0.97 \
--strand both \
--biomout ${outdir}/${sample}_otu_tax.txt \ #输出 biom 格式的结果文件
--fastapairs \ #fasta 格式文件，保存了查询序列以及对应的目标序列；
--userfields query+target+id+qcov+tcov \ #自定义输出结果的列名，从左至右分别为：查询序列 id，目标序列 id，相似度，查询序列覆盖度，目标序列覆盖度；
--userout ${outdir}/${sample}_stat.xls \ #按--userfields 定义的表头输出自定义的结果文件
以下为Usearch的使用，接去冗余那一步之后
usearch -unoise3 uniques.fa
-zotus ASVs.fa
-minsize 9
head ASVs.fa
vsearch -usearch_global filtered.fa #合并+过滤
--db ASVs.fa #合并+过滤+去冗余+聚类+去嵌合，相当于是个数据库
--id 0.99
--otutabout ASV_counts.txt
usearch -sintax ASVs.fa
-db rdp_16s_v16_sp.fa
-tabbedout ASV_tax_raw.txt
-strand both
-sintax_cutoff 0.5
网友评论