Archives

扩增子数据分析之构建OTU:UPARSE

Drive5 对扩增子测序数据分析提供了很多有益的材料,本材料由于主要介绍Usearch系列软件,很多材料就直接引自Drive5

关于如何定义"OTU", 可以直接参考: Defining and interpreting OTUs, 给出了很多有意义的信息。

"OTU" 概念来源于DNA测序之前,在数量分类学(numerical taxonomy)作为一种定量策略进行物种分类,并构建系统进化关系, 当前有更多的策略进行计算系统进化,比如16S rRNA, WGS基于核酸的方法。

在16S扩增子数据分子中,OTU被定为按照97%进行聚类,并被认为 OTU可以表述, UPARSE在构建OTU时也是采用97%这一阈值。

为什么被定义成97%,可以参考:DOTUR, a computer program for defining operational taxonomic units and estimating species richness, 其中:

    Sequences with greater than 97% identity are typically assigned to the same species,         \
those with >95% identity are typically assigned to the same genus, and those with >80% identity  \
are typically assigned to the same phylum, although these distinctions are controversial.

QIIME 1.9.x 官方网站也指出 pick_otus.py – OTU picking

An OTU threshold (default is 0.97, roughly corresponding to species-level OTUs);

UPARSE文章发表在Nature Method: UPARSE: highly accurate OTU sequences from microbial amplicon reads
在发布以来被广泛使用,截至29 Oct 2017, 已经被引用 1750 次。

这篇文章主要介绍UPARSE的算法原理和实例。

1.UPARSE 算法介绍

UPARSE 算法描述:

UPARSE的主要目的就是从已经按照丰度大小排序的序列集中获得OTU代表序列并去除嵌合体。

Uparse

上图表述了UPARSE OTU 算法, 可以解读以下几点:

1.  OTU sequence 即Centroid 序列,OTU簇代表序列,OTU序列之间 < 97% 相似度。 OTU簇成员之间≥ 97% 相似度;
2.  OTU序列在簇中丰度大于等于其它成员丰度;
3.  嵌合体序列被去除;
4.  非嵌合体序列至少要有一条序列≥ 97% 相似度,每个OTU至少包含两条序列;
5.  序列分配OTU时可能匹配 >1 个OTU, 按照序列相似度最高选择;

参考:最高丰度OTU还是最高相似度OTU? Mapping reads to OTUs

UPARSE-REF 实现逻辑:

UPARSE-REF
按照fastq-uniques结果序列按照丰度从大到小排序, 然后按照UPARSE-REF算法就是挑选OTU序列, 逻辑如下,

1. 初始化一个空的Database;
2. 序列和Database序列比对; 按照以下几种情况执行:
    a. 嵌合体序列,删除;
    b. 如空,添加到数据库;
    c. 在Database中存在序列 ≥ 97% 序列,PASS;
    d. 在Database中不存在序列 ≥ 97% 序列,添加到数据库;

从逻辑上可以看出,每条OTU序列都是丰度值从高到低选择,背后的逻辑是, 丰度高的序列最可能是真实的生物学序列,并可以用来代表具有错误(测序错误和PCR错误),嵌合体也已经被去除,因此OTU流程中不需要单独在执行去除嵌合体操作。

嵌合体序列如下图,通过比对确定模型, 并重新计算相似度,如满足≥97%,会被丢掉。

Chimers

其中UPARSE-REF关键点在于寻找最优的模型,可能是一条参考序列,也可能是多个序列组成的嵌合体模型,该方式使用了动态规划实现。

在实现过程中 singletons 序列也会被去掉,计算OTU表时大部分singletons可通过比对重新召回。

2. 命令行接口

cluster_otus 实现了UPARSE算法,主要的命令行接口如下:

-otus        : OTUs 代表序列;
-uparseout   : 描述序列如何被归类文件(分配到哪个OTU或者嵌合体)
-relabel     : OTU标识符的前缀;
-uparsealnout: UPARSE-REF模型的序列比对信息;
-minsize     : 每个OTU最小序列数, 默认2个,去除了singleton序列;

3 实例

usearch -cluster_otus uniques.fasta  -otus otus.fa -uparseout uparse.txt -relabel Otu  -uparsealnout  uparse.align.txt
usearch -otutab reads.fq -otus otus.fa -otutabout otutab.txt -mapout map.txt

整个 UPARSE 流程可以参考: OTU / denoising pipeline

获得的otus.fa需要使用 otutab 将序列比对到代表序列, 并确定OTU丰度表, 因为我们使用的未经质量控制的reads.fq序列,所以最终这些序列会分为以下5中情况:

OTU Coverage

(a) 按照97%相似度可以比对到OTUs序列的序列集合,这部分序列被认为是正确的序列,包含由于测序错误和PCR点错误造成的3%错误序列,最终被计入OTU表;
(b) 大于3%的错误的序列,未被计入OTU表;
(c) 嵌合体序列,其中一小部分会计入OTU表, 小于3%错误率会被统计如特定OTU。
(d) 具有一定错误率(大于3%)的嵌合体序列,未被技术OTU表;
(e) 真实的生物学序列,未被计入OTU表,一般为低丰度序列, 因为最小cluster_otu的 mini_size 值2或者gama值 == 4, 会导致真实的生物学序列未形成OTUs或者ZOTUs。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有。

Last Update: 2017-10-30 10:07 AM

Comments are closed.