Archives

如何使用Biostack提供的工具

我们分享的工具集都在 Github (地址) 上, 目前只提供编译后的版本, 如下图:

Biostack Github 仓库

如何下载使用工具,目前所有的工具都是在 CentOS 7编译,可能不兼容版本比较低的系统。

1. 下载工具

使用 git 直接克隆对应仓库, 比如:

git clone https://github.com/jameslz/blast_utils

如果提示没有安装 git, 请先安装git:

sudo yum install git 2. 添加环境变量

如果我们将 blast_utils 克隆到了目录: /project/tools/blast_utils, 只需要执行:

export PATH=/project/tools/blast_utils:$PATH

就可以使用这些程序了, 或者执行下面命令添加到bashrc文件,每次启动bash都会自动加载该目录:

echo export PATH=/project/tools/blast_utils:$PATH >> ~/.bashrc source ~/.bashrc

如果提示权限不够,请执行:

chmod -R 775 /project/tools/blast_utils 3. 命令执行

随意切换到任何目录,我们可以直接使用 […]

使用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序列丰度。 […]

扩增子序列样本序列丰度谱: fastx_uniques、fastx_uniques_persample、 DADA2 和 seqtk_seqtab

本文介绍扩增子序列中的一个分析需求,扩增子序列中样本序列丰度谱,即:Uniques 序列在每个样本的Tag分布, DADA2使用这个概念构建类似OTU表的特征表, QIIME2 也在使用这样特征表概念,使用这个策略,一定要有去噪 这个数据处理过程,去除测序错误造成的序列多样性。

本文只从技术层次上介绍如何产生这样的特征表,具体生物学意义可以参考下面两篇文章:

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

fastx_uniques 计算fasta/q序列中的Uniq序列,支持以下参数(详细信息见:fastx_uniques command ):

-threads :线程数; -uc :标准输出格式; -sizeout :输出簇大小信息; -relabel :重新标注序列标识符; -fastaout : fasta文件格式输出; -minuniquesize: 最小序列数/per sample 阈值,如果低于这个数值,序列不输出; -topn N : 输出按降序排列的前N条序列; -strand […]

扩增子数据分析之标准化:抽样 or 还是标准化?

扩增子测序数据会遇到样品之间数据量相差很大的情况,这是可能有数据抽平或者数据标准化概念,从技术上可认为有以下操作的方式:

1. 针对序列文件进行抽样,可以是原始测序文件,也可以是合并后的双端序列; 2. 针对OTU表进行标准化; 3. 采用DESeq2进行对原始数值进行标准化。

下面仅仅从技术操作上进行描述:

1. 序列文件进行抽样

使用seqtk samples 对原始数据进行抽样, 操作很简单:

seqtk sample -s 11 OS2_1.R1.fq.gz 5000 > OS2_1.R1.5000.fq.gz seqtk sample -s 11 OS2_1.R2.fq.gz 5000 > OS2_1.R2.5000.fq.gz

采用相同的seed值, 可以保证可重复,也可以保证R1和R2序列抽样结果可以对应。

针对合并后后序列也是同样操作:

seqtk sample -s 11 merged.fq 5000 > merged.5000.fq

这种策略虽然在概念上对数据进行重新抽样,但是会对构建OTU表构成影响,比如会人为丢掉低丰度序列, 不管uparse或者unoise对丰度都有最低需求, Uparse构建的OTU表会丢掉size<2的序列, unoise 参数默认会设置 <4个。

2. OTU表进行标准换

通常我们使用相对丰度(OTU频率表)进行描述数据,这样可以不用去考虑每个样本的测序量的大小,但是还是可以将所有的序列数(Tags数)映射到统一的测序量上,比如30000条, 可以保证相互丰度不变的情况下对数据进行标准化,该数值会对基于counts 数的计算有影响, 比如多样性指数计算等。 该操作可以通过USEARCH otutab_norm进行计算:

[…]

扩增子测序分析之多样性指数:beta 多样性

1. beta多样性介绍

Beta Diversity是对不同样品间的微生物群落构成进行比较,根据样本的OTUs丰度信息计算各种距离测度: 比如 Bray Curtis,Weighted UniFrac和Unweighted UniFrac等来评估不同样品间的微生物群落构差异。

Bray Curtis : 距离是生态学上反应群落之间差异性最常用的指标,只考虑了物种的丰度信息。 Unweighted UniFrac : 距离是基于物种(即OTUs)系统发生树进行计算的样本间的距离,只考虑了物种的有无。 Weighted UniFrac : 距离是结合OTUs的丰度信息和物种系统发生树的信息来获得的样本间的距离。

Unweighted UniFrac距离对稀有物种比较敏感,而Bray Curtis和Weighted UniFrac距离则对丰度较高的物种更加敏感,beta多样性主要反映在如何度量距离的概念。

2. 距离测度

距离的测度主要有两种形式:phylogenetic 多样性指数 和Non-phylogenetic多样性指数:

a. NON-phylogenetic 多样性指数,主要有Bray Curtis/Jaccard/euclidean等。 b. Phylogenetic 多样性指数,主要Unweighted UniFrac和weighted UniFrac

公式描述形式使用 Usearch beta diversity metric 描述。

2.1 Bray Curtis

Bray–Curtis dissimilarity, 参考wikipedia,公式很简单:

即:

Bray-Curtis = […]

扩增子数据分析之稀释曲线: alpha_div_rare 和 alpha-rarefaction

1.稀释曲线介绍

Rarefaction Curve即稀释曲线,是从样品中随机抽取一定测序量的数据,统计它们所代表物种数目(也即是OTUs数目)或多样性指数,以数据量与物种多样性来构建的曲线,以用来说明样品的测序数据量是否合理,并间接反映样品中物种的丰富程度。稀释性曲线图中,当曲线趋向平坦时,说明测序数据量渐进合理 如下图所示: 蓝色趋于饱和, 红色还一直在增加。

2.usearch alpha_div_rare

Usearch 提供了alpha_div_rare (cmd_alpha_div_rare)命令进行对OTU表进行subsampling 和计算多样性指数。 支持命令行参数:

-metric :支持的多样性指数: berger_parker、buzas_gibson、chao1、 dominance、equitability、jost、jost1. reads、richness(默认)、robbins、simpson shannon_e、shannon_2、shannon_10 -method :抽样方法,支持fast(默认), with_replacement和without_replacement。 -output :输出文件; -iters : with_replacement和without_replacement subsampling方式支持迭代数并计算平均数,并转换成最近整数。

实例:

usearch -alpha_div_rare otutab_norm.txt –method fast –metric shannon_e -output rare.txt

具体选择哪种subsample方法,结果都差不多。

USEATCH的稀释曲线和QIIME稀释曲线不同,QIIME使用counts数subsampling , Usearch使用百分比subsampling,这样的好处是所有样本都会被采样到,但是横坐标的subsampling数目不一样,会根据样本序列数目变化。

3.QIIME2 alpha-rarefaction

qiime2 数据分析流程通过 qiime diversity接口提供了分析alpha多样的各种命令: alpha-rarefaction 支持的参数:

–i-table : FeatureTable –p-max-depth: […]

扩增子数据分析之多样性指数: 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的频率; […]

扩增子测序分析之构建OTU树: USEARCH和QIIME2

扩增子数据分析中计算多样性指数可以结合系统进化树信息,比如:faith PD Tree、UniFrac 距离计算等, 本文介绍如何使用OTU序列(ZOTUs序列)构建树。

1. 距离矩阵构建树

距离在不同的场景有不同的意义:

序列相似性距离: USEARCH 定义距离 `d` = 1- identity% 衡量两条序列之间的差异, identity为两条序列比对结果中匹配的碱基与比对长度的比值。

使用 calc_distmx 计算FASTA/FASTQ距离矩阵: 支持参数:

-format :矩阵格式: tabbed_pairs(默认), square, phylip_square 或者 phylip_lower_triangle. -threads :线程数; -distmxout :矩阵输出文件; -id : 序列最低相似度;

使用 cluster_aggd可对输出距离矩阵构建树,构建树的方式可以通过linkage指定,支持参数:

-treeout : 输出构建的树(Newick格式); -linkage : 等级聚类模型,min|max|avg

实例:

usearch -calc_distmx otus.fa -distmxout dist.txt -format square usearch -cluster_aggd mx.txt -treeout […]

扩增子数据分析之去噪: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%的非冗余序列描述种或者株,首先要解决的问题就是鉴定哪些序列簇来自同一个生物学序列模板,这个就是我们说的去噪。

去噪 […]