Archives

扩增子序列分析之物种分类: 注释到种水平,也许可以!

一、引言

DADA2 升级至 1.6 版本, 并正式支持 Species assignment, 使用了精确匹配模式判别是否可以分类到种水平,当然种以上分类水平还是要使用RDP算法或者Sintax算法。

DADA2 Species assignment 逻辑很简单:

预处理16S库文件(去除未分类的序列或者分类不明确的序列, 保留属和种水平分类信息), DADA2 提供了rdp_species_assignment_16.fa.gz和silva_species_assignment_v128.fa.gz 两种库。 序列精确比对,DADA2使用了ShortRead进行序列比对; 判断,如果所有有效匹配(可设正向和反向选项)在 属+种 分类都一致,那么是可以直接分配该种分类等级, DADA2 提供了允许分配多个种的选项,不推荐使用这个模式。

代码段:

assignSpecies <- function(seqs, refFasta, allowMultiple=FALSE, tryRC=FALSE, verbose=FALSE) { # Define number of multiple species to return if(is.logical(allowMultiple)) { if(allowMultiple) keep <- Inf else keep <- 1 } else { […]

使用Usearch进行扩增子数据分析:The Full Story

准备工作:

添加辅助程序路径至环境变量:

export PATH=$PATH:/biostack/tools/pipelines/div_seq-0.2.3/bin/utils export PATH=$PATH:/biostack/labs/seqtk_utils/release export PATH=$PATH:/biostack/labs/tabtk_utils/release export PATH=$PATH:/biostack/tools/microbiome/KronaTools-2.7/bin

准备好原始数据,样品信息表、以及引物信息。 软件需要安装:R、seqtk、usearch、seqtk_utils、tabtk_utils、biom以及fasttree、 Trimal, mafft、KronaTools以及phylommand

1. 数据整体质量评估

参考: 获取fasta/q文件的碱基组成和长度信息: seqtk_info 和 seqtk_fqchk

mkdir quality seqtk_fqchk data/raw.R1.fastq.gz >quality/raw_R1.tsv fqchk_base.R quality/raw_R1.tsv quality/raw_base_R1.pdf R1 fqchk_qual.R quality/raw_R1.tsv quality/raw_qual_R1.pdf R1

可以生成质量和碱基分布图:

2. 数据拆分

参考: 快速对扩增子序列进行数据拆分: seqtk_demultiplex

mkdir demultiplex cut -f1,2 mapping_file.txt | sed ‘s/,/\t/’ >demultiplex/barcode.txt seqtk_demultiplex -b demultiplex/barcode.txt -1 […]

扩增子测序分析之去噪 – DADA2及其平台化

引言

DADA2(Divisive Amplicon Denoising Algorithm)的核心 去噪 然后使用 Amplicon Sequence Variants (ASVs)概念构建类OTU表;

具体算法可以参考下面几篇文章:

Denoising PCR-amplified metagenome data DADA2: High-resolution sample inference from Illumina amplicon data Bioconductor Workflow for Microbiome Data Analysis: from raw reads to community analyses

二、核心算法

DADA2 的核心算法:

1. 错误率模型计算,采用selfConsist无监督学习模型; 2. Rate λ_ji: 通过错误率模型,衡量扩增子序列i是否来在来自j; 3. 丰度P-value值,使用丰度信息衡量:该丰度的扩增子序列是否可以通过错误模型解释,使用了一个经验P值来判断,比如 OMEGA_A = 1e-40 4. 分裂算法( divisive partitioning)形成扩增子簇,簇心定义为denoise序列,簇大小定义denoise序列丰度。 […]

扩增子数据分析之多样性指数: alpha多样性

多样性指数(Diversity index)和计算公式可以见: wikipedia

Alpha多样性(Alpha Diversity)是对某个样品中物种多样性的分析,包含样品中的物种类别的多样性——丰富度(Richness)和物种组成多少的整体分布——均匀度(Evenness)两个因素,通常用Richness,Chao1,Shannon,Simpson,Dominance和Equitability等指数来评估样本的物种多样性。

丰富度指数

Richness, Chao1,Shannon三个指数是常用的评估丰富度的指标,数值越高表明样品包含的物种丰富度就越高。

Richness指数: 指样本中被检测到的OTU量; Chao1指数 : 通过低丰度OTUs来进一步预测样品中的OTUs数量; Shannon指数 : 计算考虑到样品中的OTUs及其相对丰度信息, 通过对数(如以2为底的shannon_2,以自然对数为底的shannon_e 以10为底的shannon_10)转换来预测样品中的分类多样性。

均匀度指数

Simpson,Dominance和Equitability三个指数是常用的评估均匀度的指标。

Simpson指数 : 表示随机选取两条序列属于同一个分类(如OTUs)的概率(故数值在0~1之间), 数值越接近1表示表明OTUs的丰度分部越不均匀; Dominancez指数 : 取值为1-Simpson,表示随机选取两条序列属于不同分类(如OTUs)的概率; Equitability指数: 根据Shannon指数值计算,当其值为1时表明样品中的物种丰度分布绝对均匀, 而其值越小这表明物种丰度分布呈现出越高的偏向。

汇总表:

指数 单位 计算方式 richness OTUs 样本中至少包含一条序列的OTU数目 chao1 OTUs N + S^2 / (2D^2),其中N为OTU个数, S为丰度为1的OTUs个数,D为丰度为2的OTUs数目; shannon_2 bits sum(f), 对所有OTU频率计算p*log(p,2)和, p为OTU的频率; shannon_e nats sum(f), 对所有OTU频率计算p*log(p,e)和, p为OTU的频率; […]

扩增子数据分析之去噪:UNOISE

Error Cloud 概念模型:

子图c刻画了 Error Cloud 概念,每个圆圈代表一个Unique序列,圈大小描述序列数,绿色圆圈表述真实生物学序列, 黄色圈描述了替换突变导致的变体,中间黑线为编辑距离(两条序列差异的个数, 大部分为d=1)。. 该图来自:Interpreting 16S metagenomic data without clustering to achieve sub-OTU resolution

先祭出两篇文献:

The 97% threshold is much too low, especially for V4 Exact sequence variants should replace operational taxonomic units in marker-gene data analysis

当前开始认为使用 100% 的相似度去度量生物学种或者株的更合适, 97%是一个在数据量不是很充足的条件下产生的结论,Usreach的UNOISE和DADA2都采用了这样的概念。

既然使用100%那问题来了, 先前97%通过聚类可以掩盖一些噪音, 这些噪音包括:

测序错误; PCR错误(点错误和嵌合体); 菌株之间的变异;

现在使用100%的非冗余序列描述种或者株,首先要解决的问题就是鉴定哪些序列簇来自同一个生物学序列模板,这个就是我们说的去噪。

去噪 […]

扩增子数据分析之构建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 […]

扩增子数据分析之聚类:UCLUST

UCLUST 使用USEARCH作为序列比对引擎并对序列进行聚类, 聚类采用了一种贪婪算法(是一种在每一步选择中都采取在当前状态下最好或最优(即最有利)的选择,从而希望导致结果是最好或最优的算法 wikipedia)。

描述 UCLUST的文章发表在Bioinformatics杂志: Search and clustering orders of magnitude faster than BLAST

QIIME 1.9.x 版本的构建OTU使用的就是UCLUST算法, 本文主要介绍UCLUST算法以及想对于应等级聚类概念。

1. UCLUST算法介绍

UCLUST 可对氨基酸序列和核酸序列进行聚类,其算法原理如下:

UCLUST对序列进行聚类,寻找簇集,尝试满足下面个条件, T为指定相似度阈值:

(1) 所有的簇心之间相似度 < T (2) 簇成员之间相似度 ≥ T

UCLUST的簇心选择和输入文件的序列顺序有很大关联,序列可以按照序列长度进行排序,也可以按照序列的丰度进行排序,设计的逻辑是序列在文件中先出现先获得作为簇心的机会,并且每个簇只有一条代表序列。

实现逻辑:

(1). 初始化空的database; (2). 将序列文件中序列按照顺序开始放入database, 根据以有几种情况进行选择: a. Database 为空, 放入database; b. 和Databases 序列比对,判断相似度阈值: 如果相似度 ≥T 置入对应的簇。 如果相似度 <T 添加到database, 该序列设置为簇心; […]

去除嵌合体: uchime_ref和uchime3_denovo

嵌合体从序列组成来说是指一条序列由两条或者更多条生物学序列连接而成(见下面示意图),扩增子测序实验中,主要由于PCR阶段不完全延伸形成,形成的嵌合体会在随后PCR过程被放大。具体解释可参考 chimera_formation。

嵌合体示意图:

下面讨论如何去除这些嵌合体,以及什么场景采用什么策略。

USEARCH系统去除嵌合体策略主要有:

使用嵌合体参考库 参考:uchime2_ref 从头去除嵌合体 参考:uchime3_denovo OTU构建集成: 参考: uparse 和 unoise3

请参考: UCHIME2 算法: UCHIME2: Improved chimera detection for amplicon sequences

1. 使用嵌合体参考库

命令行接口:

usearch -uchime_ref reads.fasta -db silva.udb -uchimeout out.txt -strand plus -threads 20 -mode balanced

主要参数解释:

-db: 参考库, 推荐SILVA; -strand: 目前只支持 plus,所以序列最好是调整好方向的; -uchimeout: 嵌合体输出文件; -threads: 线程数; -mode: 可以使用的模型,支持 high_confidence、specific、balanced、sensitive等几种模式。\ […]