16S扩增子之Vsearch使用

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

4. 去冗余与聚类生成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

整理自https://mp.weixin.qq.com/s/MTZFWIqr1dqO-fwUQ8Pb4g