病原微生物扩增子数据分析实战（三）：vsearch软件鉴定物种组成

病原微生物基因检测的两大核心任务是物种组成和功能组成的鉴定，而扩增子测序的首要目的是找到致病的细菌或者病毒，即鉴定物种组成。

经过质量控制，测序数据中已经不包含非生物的异常序列了，下面我们用vsearch[1]软件完全后续分析。

1.去重（dereplication）

同一对引物的扩增产物，理论上应该是完全一样的，这些冗余的信息会造成比较大的运算负担，因此需要先去冗余，相同的序列只保留一条就好了。

vsearch 去冗余有以下三种模式：

• 全长模式：两条序列长度完全相同，碱基完全一样，总之，序列要完全一模一样才认为是重复
• 一致模式：在全长模式的基础上，进一步要求序列的名称要完全相同
• 前缀模式：只要一条序列是另一条序列的前缀，也就是一条序列与另一条序列的前面部分相同，就可以认为它们是重复。

可以看出，这三种模式的严格程度依次是：一致模式 > 全长模式 > 前缀模式。

我们用全长模式去除冗余：

vsearch --derep_fulllength ${outdir}/${sample}_merge.fastq.gz --output ${outdir}/${sample}_derep.fa --minuniquesize 8 --minseqlength 100 --strand both --sizeout --fasta_width 0

--derep_fulllength，后面跟通过质控的 fastq 文件；

--output，输出文件，fasta 格式；

--minuniquesize，最低丰度值，低于该丰度的序列会被过滤掉；

--minseqlength，最低长度值，低于该长度的序列会被过滤掉；

--strand，当判断两条序列是否一致时，默认只考虑正链plus，both表示考虑正反两个方向都考虑；

--sizeout，在结果文件中序列名称后面添加丰度信息；

--fasta_width，限定 fasta 结果文件中每条序列在一行中最多显示的字符数，默认是 80，0 表示不做限制；

2.降噪（denoise）

按 97%相似度对序列进行聚类曾经是扩增子序列分析的金标准，但这有一个问题，就是物种只能鉴定到属或种，因为有些物种之间的相似度接近 99%，甚至有的菌株之间只有 SNP 的差别。因此现在流行的做法是进行降噪处理，通俗地说就是按 100%相似度进行聚类，从而使鉴定达到种或株的水平。

vsearch --cluster_unoise ${outdir}/${sample}_derep.fa --centroids ${outdir}/${sample}_cent.fa --consout ${outdir}/${sample}_cons.fa --minsize 8 --strand both --sizeout --fasta_width 0 --clusterout_sort

--cluster_unoise，上一步去重后的 fasta 文件；

--centroids，fasta 结果文件，包含每一个聚类中的种子序列；

--consout，fasta 结果文件，包含每一个聚类的一致性序列；

--minsize，降噪最低的序列丰度要求，默认是 8 条；

--strand，当判断两条序列是否一致时，默认只考虑正链plus，both表示考虑正反两个方向都考虑；

--sizeout，在结果文件中序列名称后面添加丰度信息；

--fasta_width，限定 fasta 结果文件中每条序列在一行中最多显示的字符数，默认是 80，0 表示不做限制；

--clusterout_sort，结果文件中序列的顺序默认是按其在输入文件中的顺序，设定该参数则是按照降噪后序列的丰度排序；

3.去嵌合体

测序数据中可能存在嵌合体序列，每个嵌合体至少来源于两个或两个以上的扩增子模板，因而需要事先去除。

去嵌合体有两种方法：一是基于参考数据库，一是基于无参考数据库的 denovo 方法。

vsearch 基于 denovo 方法去除嵌合体的命令：

vsearch --uchime3_denovo ${outdir}/${sample}_cons.fa --nonchimeras ${outdir}/${sample}_nochimeras.fa --chimeras ${outdir}/${sample}_chimeras.fa --uchimealns ${outdir}/${sample}_chimeras.aln --sizeout --fasta_width 0

--uchime3_denovo，上一步去重后获得的一致性序列文件；

--nonchimeras，不含有嵌合体的结果文件；

--chimeras，含有嵌合体的结果文件；

--uchimealns，以人类易于阅读的形式呈现嵌合体与其两个亲本进行比对的结果文件；

--sizeout，在结果文件中序列名称后面添加丰度信息；

--fasta_width，限定 fasta 结果文件中每条序列在一行中最多显示的字符数，默认是 80，0 表示不做限制；

4.创建 OTU 表

OTU（operational taxonomic units），即可操作分类单元，通过 vsearch 创建 OTU 表，可以看到每个 OTU 的丰度情况。

a.拷贝 OTU 文件：经过去重、降噪、去嵌合体之后，获得的无嵌合体的 fasta 文件可作为最终 OTU 结果

cp ${outdir}/${sample}_nochimeras.fa ${outdir}/${sample}_otu.fa

b.格式转换：将 fastq 格式的质控文件转换成 fasta 格式，以下命令只做格式转换，没有进行任何过滤操作

vsearch --fastx_filter ${outdir}/${sample}_merge_clean_final.fastq.gz --fastaout ${outdir}/${sample}_merge_clean_final.fa

c.创建 OTU 表：本质上是将经过质控的序列比对到 OTU 序列上，从而可以得出每一个 OTU 序列的丰度

vsearch --usearch_global ${outdir}/${sample}_merge_clean.fa --db ${outdir}/${sample}_otu.fa --id 0.97 --query_cov 0.97 --strand both --otutabout ${outdir}/${sample}_otutab.txt --samheader --samout ${outdir}/${sample}.sam

--usearch_global，后面跟经过质控的 fastq 文件，fasta 格式；

--db，OTU 文件，fasta 格式；

--id，相似度阈值：当查询序列与目标序列之间的相似度达到多少时，才算比对上；

--query_cov，覆盖度：满足相似度的情况下，同时要求查询序列的覆盖度达到多少；

--strand，当判断两条序列是否一致时，默认只考虑正链plus，both表示考虑正反两个方向都考虑；

--outtabout，输出 OTU 表。

--samheader，当要保存 SAM 文件时，默认不包含 SAM 文件头，设定此参数即会包含；

--samout，输出比对的 SAM 格式结果文件，可用于后续 IGV 软件可视化查看比对情况。

5.物种鉴定

有了样本的代表性序列后，想知道样本的物种组成，可将这些序列比对到参考数据库，通过相似性及覆盖度进行判断。

vsearch --usearch_global ${outdir}/${sample}_otu.fa -db $DB/amplicon.fa --id 0.97 --query_cov 0.97 --strand both --biomout ${outdir}/${sample}_otu_tax.txt --alnout ${outdir}/${sample}.aln --blast6out ${outdir}/${sample}_blast.xls --fastapairs ${outdir}/${sample}_pairs.fa --notmatched ${outdir}/${sample}_notmatched.fa --userfields query+target+id+qcov+tcov --userout ${outdir}/${sample}_stat.xls

--usearch_global，OTU 文件，fasta 格式；

--db，参考序列库，fasta 格式；

--id，相似度阈值：当查询序列与目标序列之间的相似度达到多少时，才算比对上；

--query_cov，覆盖度：满足相似度的情况下，同时要求查询序列的覆盖度达到多少；

--strand，当判断两条序列是否一致时，默认只考虑正链plus，both表示考虑正反两个方向都考虑；

--biomout ，输出 biom 格式的结果文件；

--alnout，输出比对结果文件；

--blast6out，输出 blast 6 格式的比对结果文件；

--fastapairs，fasta 格式文件，保存了查询序列以及对应的目标序列；

--notmatched，没有搜索到目标序列的查询序列集合；

--userfields，自定义输出结果的列名，从左至右分别为：查询序列 id，目标序列 id，相似度，查询序列覆盖度，目标序列覆盖度；

--userout，按--userfields 定义的表头输出自定义的结果文件。

上述命令得到一个自定义的结果文件{outdir}/{sample}_stat.xls，根据其中的相似度、覆盖度以及前面的 OTU 表{outdir}/{sample}_otutab.txt，就可以知到物种的组成情况，以及它们的丰度信息。

现在，通过 3 篇文章，一套病原微生物扩增子物种鉴定的流程就建立起来了。

「往期精彩」

病原微生物扩增子数据分析实战（二）：fastp 软件进行质量控制

病原微生物扩增子数据分析实战（一）：bcl2fastq 软件完成数据拆分

Reference

[1]

vsearch: https://github.com/torognes/vsearch