Archives

16S扩增子组成谱预测功能谱:本地快速执行otutable_norm和otutable_melt以及otutable_kann

当前可以使用16S扩增子组成谱预测功能谱的工具主要有:PICRUStTax4funPiphilin, 本文主要使用PICRUSt的预测文件,进行重新实现。

PICRUSt 采用了一种ASR(ancestral state reconstruction, 可参考1985年的一篇文章:Phylogenies and the Comparative Method)方法重构物种的16S拷贝数和基因组的基因功能信息(使用KEGG、COG描述),基于系统进化的逻辑是同一个clade分支下的基因组(物种)具有相似的形状,这样就可以使用已测序基因组预测未测序基因组的功能空间,主要依赖参数为系统进化树分支距离。

PICRUSt

图来自PICRUSt文章

原理图分拆解读:

1. 使用Greengene 16S序列构建系统进化树,去除未测序基因组;
2. 预测所有祖先状态的性状(16S拷贝数、基因功能);
3. 预测未知基因组节点的性状(16S拷贝数、基因功能);
4. 标准化,使用预测的拷贝数对OTU丰度进行标准化;
5. 基因功能空间预测(KO、COG)丰度表;

其实关键算法在1、2、3部分。 这篇文章我们直接使用PICRUSt根据 Greengene 13.5 进行的预测结果进行一些数据分析。

需要的文件为双端合并后的序列数据(Usearch sample=’xxx’ 注释文件),我们会在下面一步一步演示如何操作。 准备工作:

下载PICRUSt预计算拷贝数和KO表: 16S_13_5_precalculated.tab.gzko_13_5_precalculated.tab.gz 下载Greengene 97OTU代表序列: 97_otus.fasta.gz

如何使用这些程序,请移步: 如何使用biostack 提供的命令行工具

1. 构建OTU表

虽然有文章讨论过使用closed- 和 open-reference 构建 OTUs 有多不准确,如下:

Accuracy of microbial community diversity estimated by closed- and open-reference OTUs.
Edgar RC.
PeerJ. 2017 Oct 4;5:e3889. doi: 10.7717/peerj.3889. eCollection 2017.
PMID: 29018622

如果有兴趣看看问题在哪里, 可以移步:closed_ref_problems, 这不是这篇文章的主题,因为我们使用基于Greengene 13.5 的PICRUSt预先计算的结果,所以只能选择使用Greengene作为参考构建OTU表。

这里使用Usearchclosed_ref 构建closed OTU表;

closed_ref 命令行接口:

usearch -closed_ref reads.fastq -db gg97.fa -otutabout otutab.txt -strand plus -tabbedout closed.txt

下面简单解读:

reads.fastq: 为具有 'sample=xxx' 注释的16S扩增子序列集合。
-db        : 为greengene 97%聚类 OTU代表序列。
-strand    : plus或者both,搜索查询序列的条件,只都搜索正链还是正负链都搜索。
-otutabout : OTU表输出;
-tabbedout : closed_ref 比对report文件;
-id        : 最低相似度要求,默认参数97%;

下载Greengene 97% OTU代表序列文件:

wget https://github.com/biocore/qiime-default-reference/blob/master/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta.gz

执行构建OTU表:

gunzip 97_otus.fasta.gz
usearch -closed_ref reads.fastq -db 97_otus.fasta -otutabout otutab.txt \
        -strand plus -tabbedout closed.txt -threads 20

2. 16S OTU表 标准化

上一步我们已经获得了OTU表信息,下一步就是使用Greengene的16S 拷贝数进行标准化OTU表; 我们使用 `otutable_norm` 程序实现标准化工作

命令行接口:

$ otutable_norm
Usage: otutable_norm <copy_number> <otutab>
version: 0.0.1

otutable_norm 接受两个参数: 拷贝数表和OTU表;

执行:

wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust_precalculated_v1.1.2/13_5/16S_13_5_precalculated.tab.gz
otutable_norm  16S_13_5_precalculated.tab.gz  otutab.txt >otutab_norm.txt

3. 合成功能丰度表

这一步使用 otutable_melt 完成合成功能谱的工作。 命令行接口:

$otutable_melt
Usage: otutable_melt [options] <feature-matrix> <otu_table>
Options:
  -r  round value to nearest intege;
  -v  print version number

otutable_melt 接受两个参数:PICRUSt兼容的特征表和OTU表(QIIME兼容的OTU表),这里使用上一步产生标准化后的OTU表, otutable_melt也接受 Bugbasedefault_traits_precalculated.txt数据;

-r: 四舍五入,转换成整型输出。

执行:

wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust_precalculated_v1.1.2/13_5/ko_13_5_precalculated.tab.gz
otutable_melt   ko_13_5_precalculated.tab.gz otutab_norm.txt  >ko_profiling.txt

或者:

otutable_melt   <(pigz -dk -c  -p 20  ko_13_5_precalculated.tab.gz)  otutab_norm.txt  >ko_profiling.txt

现在我们已经可以获取KEGG的入门券了,后面我们可以使用KO关联KEGG的很多重要信息, 比如 Pathway,Module,Enzyme等。

4. otutable_kann

KO(KEGG Ortholog)是KEGG唯一一个通过序列和功能关联的中间点,我们获得KO的丰度信息, 可以通过合并同一个Pathway或Module的KO丰度,计算Pathway和Module的丰度。

实现该功能我们使用 otutable_kann实现;

命令行接口:

$ otutable_kann
Usage: otutable_kann [options] <feature-map> <otu_table>
Options:
  -r  round value to nearest intege;
  -v  print version number

otutable_kann接受两个参数, KEGG的KO映射文件和OTU表文件(KEGG或COG的丰度表这里统一为otu_table表)。 文件为两列:

第一列: KO Number 第二列: KO 映射表;

该文件可以通过KEGG的API服务下载,比如:

curl http://rest.kegg.jp/link/pathway/ko | \
sed 's/ko://; s/path://' | grep "map"  >ko_pathway.tab
curl http://rest.kegg.jp/link/module/ko |   \
sed 's/ko://; s/md://'  >ko_module.tab

下载注释文件:

curl http://rest.kegg.jp/list/module | sed 's/md://' >module-definition
curl http://rest.kegg.jp/list/pathway | sed 's/path://' >pathway-definition
curl http://rest.kegg.jp/list/ko | sed 's/path://' >ko-definition

执行(以Module为例子):

otutable_kann   ko_module.tab  ko_profiling.txt >module_profiling.txt

如果需要注释信息需要使用 tabtk_definition 添加:

otutable_kann   ko_module.tab  ko_profiling.txt  | \
tabtk_definition -d " " module-definition  - >module_profiling.txt

5. 总结

到这里我们可以很简单高效完成共物种组成谱到功能谱的预测,我们提供的otutable_normotutable_meltotutable_kann 都比较通用,组合起来可以解决很多问题。

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

Last update:2017-11-14 7:46 PM

如何使用Biostack提供的工具

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

ScreenClip.png-13.6kB

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. 命令执行

随意切换到任何目录,我们可以直接使用 blast_utils 目录下应用程序,比如:

$ cd

$ blast_hsp
Usage: blast_hits  [options]  <blast|->
Options:
  -b DOUBLE  MIN bit score, default: [60]
  -e DOUBLE  MAX E-value, default: [0.001]
  -i DOUBLE  MIN identity, default: [0]
  -v print version number

4. 其它问题

如果可执行文件在仓库的 bin 目录, 添加目录时请包括bin目录, 比如:

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

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

Last update:2017-11-15 3:02 PM

fasta/q序列抽取和排序: seqtk_reoder

一、使用场景

给定fasta/q文件和一个序列列表,按照给定的序列顺序抽取序列,如不存在则忽略,先前的 seqtk_subseq( 如何根据列表快速提取fasta/q序列子集/补集:seqtk_subseq ) 抽取出序列按照标识符在序列文件中出现的顺序输出,有时我们需要按照给定的标识符输出序列。

seqtk_reoder使用reoder作为名字输出就是指序列输出顺序按照指定顺序输出。

二、工具介绍

通过seqtk_info 查看测试序列信息:

$seqtk_info  shear.fastq
#sequence       base    min_len max_len avg_len
47525   23905075        503     503     503.00

一共47525序列;

2.1 seqtk_split

命令行接口:

$ seqtk_reorder
Usage: seqtk_reorder <fasta/q> <list>
version: 0.0.1

seqtk_reorder 接受两个参数, 序列文件列表文件, 支持stream流模式;

实例

$cat  shear.idx | seqtk_reorder  shear.fasta -

$cat  shear.idx | seqtk_reorder  shear.fastq.gz  -

$cat  shear.idx | seqtk_reorder  <(gunzip  -c shear.fastq.gz)  -

该程序设计使用hash map实现存储目标序列,所以不支持数据库比较大的序列抽取, seqtk_subseq 只存储列表,所以适用范围更广。

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

Last update:2017-11-14 4:08 PM

fasta/q序列拆分三剑客: seqtk_split、seqtk_partition 和 seqtk_shred

一、使用场景

本部分内容主要解决以下几个问题:

  1. 对指定fasta/q序列文件拆分成每份N条序列问题,比如执行InterProScan。
  2. 对指定fasta/q序列文件拆分成指定N份问题,比如我们有40核心,程序不支持多线程(或者支持的不理想),那拆分40份多进程执行, 比如HMMER
  3. 对指定fasta/q序列文件序列拆分,每条序列保留最大l长度,并保证分割非序列存在N个碱基交叠。使用的场景多,比如 454 Newbler拼装软件,不支持很长度序列,如果我们对不同拼装版本的组装结果混合组装,将序列撕碎(保留比较大的overlap)再进行组装,就没有问题了。

二、工具介绍

以下三个程序分别解决以上三个问题。

通过seqtk_info 查看测试序列信息:

$seqtk_info  shear.fastq
#sequence       base    min_len max_len avg_len
47525   23905075        503     503     503.00

一共47525序列;

2.1 seqtk_split

命令行接口:

$ seqtk_split
Usage: seqtk_split [options] <fasta/q> <size> <prefix>
Options:
  -t  fastx file suffix: default: [fasta];
  -v  print version number

seqtk_split 接受三个参数, 序列文件序列大小(条数)输出文件前缀
可以通过 -t 指定输入文件类型, 比如 fasta或者fastq

实例

$mkdir split

$seqtk_split -t fastq  shear.fastq   10000   split/split

$ls split
split_1.fastq  split_2.fastq  split_3.fastq  split_4.fastq  split_5.fastq

seqtk_split会返回分拆的文件个数,方便自动化批量处理。

2.2 seqtk_partition

命令行接口:

$ seqtk_partition
Usage: seqtk_partition [options] <fasta/q> <partition: 10> <prefix>
Options:
  -t  fastx file suffix: default: [fasta];
  -v  print version number

seqtk_partition 接受三个参数, 序列文件分割数输出文件前缀
可以通过 -t 指定输入文件类型, 比如 fasta或者fastq

实例

$mkdir partition

$seqtk_partition -t fastq  shear.fastq   3   partition/partition

$ls partition
partition_1.fastq  partition_2.fastq  partition_3.fastq
2.3 seqtk_shred
$ seqtk_shred
Usage: seqtk_segment <fasta> <length> <overlap>

seqtk_shred 只支持fasta文件, 接受三个参数, 序列文件分割成序列长度序列交叠长度

实例

$seqtk seq -A  shear.fastq  >shear.fasta

$seqtk_shred shear.fasta 300 200 >shear.shred.fasta

输出序列会根据位置添加_1, _2, _3…

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

Last update:2017-11-14 3:11 PM

使用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

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

raw_qual_R1.png-298.2kB
raw_base_R1.png-711.2kB

2. 数据拆分

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

mkdir demultiplex
cut -f1,2 mapping_file.txt | sed  's/,/\t/'  >demultiplex/barcode.txt
seqtk_demultiplex -b  demultiplex/barcode.txt  -1  data/raw.R1.fastq.gz   -2  data/raw.R2.fastq.gz  -d demultiplex -l 5

3. 双端序列合并

参考: 双端序列合并:fastq_mergepairs

单个样本可以这样:

mkdir mergepairs
usearch -fastq_mergepairs demultiplex/D1_R1.fastq \
        -reverse demultiplex/D1_R2.fastq      \
        -fastqout mergepairs/D1.fastq         \
        -fastqout_notmerged_fwd mergepairs/D1.notmerged_fwd.fastq \
        -fastqout_notmerged_rev  mergepairs/D1.notmerged_rev.fastq \
        -log  mergepairs/D1.log \
        -threads 24 \
        -fastq_minmergelen 100 \
        -fastq_maxmergelen 500 \
        -fastq_maxdiffs 10  \
        -fastq_pctid 80 \
        -fastq_trunctail 2

可使用Bash进行循环迭代:

cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{usearch -fastq_mergepairs demultiplex/$F[0]\_R1.fastq -reverse demultiplex/$F[0]\_R2.fastq -fastqout mergepairs/$F[0].fastq -fastqout_notmerged_fwd mergepairs/$F[0].notmerged_fwd.fastq  -fastqout_notmerged_rev  mergepairs/$F[0].notmerged_rev.fastq -log  mergepairs/$F[0].log -threads 24 -fastq_minmergelen 100 -fastq_maxmergelen 500  -fastq_maxdiffs 10 -fastq_pctid 80 -fastq_trunctail 2 \n}' |bash

4. 引物识别和切除引物

参考: 对扩增子序列进行引物匹配: usearch 和 primersearch 切除引物序列获取扩增目标序列:seqtk_pcr_strip

mkdir primer_match pcr_strip
usearch -search_pcr mergepairs/D1.fastq \
        -db  primers.txt -strand both   \
        -threads 5 -maxdiffs 2          \
        --pcrout primer_match/D1.txt -log primer_match/D1.log ;

seqtk_pcr_strip -l 4  primer_match/D1.txt  mergepairs/D1.fastq  >pcr_strip/D1.fastq

通过bash可以批量执行:

cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{usearch -search_pcr mergepairs/$F[0].fastq -db  primers.txt -strand both -threads 5 -maxdiffs 2  --pcrout primer_match/$F[0].txt -log primer_match/$F[0].log ; \n}' |bash

cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{seqtk_pcr_strip -l 4 primer_match/$F[0].txt  mergepairs/$F[0].fastq  >pcr_strip/$F[0].fastq ; \n}' |bash

5. 添加序列注释信息和pool样本

参考: 添加和去除 fasta/q 序列标识符注释信息: seqtk_labels和fastx_strip_annots

mkdir pools
cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{seqtk rename  pcr_strip/$F[0].fastq  $F[0]_ | seqtk_labels -  sample=$F[0]  >> pools/seq.fq ; \n}' |bash

此时序列标识符为:

@D1_1;sample=D1;

6. 质量控制

参考: NGS数据质量过滤:Usearch fastq-filter、seqtk trimfq和seqtk_filter

usearch -fastq_filter pools/seq.fq   \
        -fastq_maxee 1.0    \
        -fastaout  pools/filtered.fa -relabel Filt

7. fastx_uniq 去冗余

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

usearch -fastx_uniques  pools/filtered.fa \
        -relabel Uniq                 \
        -sizeout -fastaout  pools/derep.fa \
        -uc pools/derep.uc      \
        -log pools/derep.log; 

8. UPARSE

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

mkdir uparse
usearch -cluster_otus pools/derep.fa -otus uparse/otus.fa -relabel OTU_
usearch -otutab  pools/seq.fq  -otus uparse/otus.fa -otutabout uparse/otutab.txt

9. Unoise

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

mkdir unoise
usearch -unoise3 pools/derep.fa -zotus unoise/otus.fa  -minsize 4

10. 物种分类 sintax

参考: 使用SINTAX对16S序列进行分类:SINTAX 使用 rdp_16s_v16_sp 进行预测分类水平。

mkdir sintax
usearch -sintax uparse/otus.fa \
    -db  /biostack/databases/sintax/rdp_16s_v16_sp.udb  \
    -threads 24  -strand both  \
    -sintax_cutoff 0.8  -tabbedout sintax/otu.sintax

11. OTU注释

参考: 计算OTU表相对丰度:otutable_freq
对OTU表进行注释:otutable_annotation

mkdir otus
tabtk_row_reorder uparse/otutab.txt  <(seqtk_names uparse/otus.fa) > otus/otutab.txt
cut -f1,4  sintax/otu.sintax |  otutable_annotation  -  otus/otutab.txt  >otus/otutab_ann.txt
otutable_freq otus/otutab.txt  >otus/otutab_freq.txt
cut -f1,4  sintax/otu.sintax |  otutable_annotation  - otus/otutab_freq.txt >otus/otutab_freq_ann.txt
tabtk_xlsx  otus/otus.xlsx  otutab:otus/otutab.txt  \
                            otutab_freq:otus/otutab_freq.txt \
                            otutab_ann:otus/otutab_ann.txt  \
                            otutab_freq_ann:otus/otutab_freq.txt

biom convert  -i otus/otutab_ann.txt  -o  otus/otutab.biom  \
        --table-type "OTU table"  \
        --to-json  --process-obs-metadata taxonomy

12. 物种分类 krona 可视化

参考: OTU表生成 Krona 输入文件:otutable_krona

mkdir krona
otutable_krona  otus/otutab_ann.txt D1 >krona/D1.txt
ktImportText   -o krona/kron.html  krona/D1.txt,D1

现在应该可以看到:

image.png-493.8kB

13. 物种分类可视化

参考: 生成不同水平系统分类丰度表:otutable_abundance 和 otutable_level

mkdir taxonomy
otutable_level  -l g  otus/otutab_freq_ann.txt  >taxonomy/genus.txt

这是可以获得每个样本在不同分类水平的丰度情况,比如在属水平。 选择前30进行数据可视化:

taxon_table_rank taxonomy/genus.txt  30  >taxonomy/genus.30.txt
barplot.R  taxonomy/genus.30.txt   taxonomy/genus.30.pdf  genus.30
pdf2png taxonomy/genus.30.pdf

genus.30.png-92.7kB

也可以绘制热图。

14. 构建系统进化树

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

mkdir phylogeny
mafft  --auto  uparse/otus.fa   >phylogeny/otus.aln
trimal  -in phylogeny/otus.aln  -out  phylogeny/otus.trimal.fasta  -automated1
fasttree  -nt  phylogeny/otus.trimal.fasta >phylogeny/otus.nwk
treebender  --clear_internal_node_labels phylogeny/otus.nwk | treebender  --mid_point_root  >phylogeny/rooted.nwk

15. alpha多样性

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

mkdir alpha
usearch -alpha_div uparse/otutab.txt -output alpha/alpha.txt
alpha_table  alpha/alpha.txt  mapping_file.txt alpha/diversity
estimators_stats.R  alpha/diversity/richness.txt  alpha/diversity/richness_boxplot.pdf richness

可以绘制多样系指数的boxplot, 如下:

richness_boxplot.png-18kB

16. 稀释曲线

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

mkdir rare_curve
usearch -alpha_div_rare  otus/otutab.txt   --method fast --metric  shannon_e  -output  rare_curve/rare.txt
rarefaction_curve.R  rare_curve/rare.txt  mapping_file.txt   rare_curve/rare.pdf

可绘制如下图:

rare.png-73.5kB

17. beta多样性

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

mkdir beta
usearch -beta_div uparse/otutab.txt -metrics bray_curtis,jaccard -filename_prefix /project/amplicon/beta/distmatrix/
usearch -beta_div uparse/otutab.txt -tree phylogeny/rooted.nwk -metrics unifrac,unifrac_binary -filename_prefix /project/amplicon/beta/distmatrix/

到此,我们使用USEARCH 完成了大部分的扩增子生物信息数据分析工作, 后续主要是多元统计以及可视化分析和功能预测, 比如,可视化组内的距离和组间的距离是否有差异(比如使用wilcox.test, 可以产生下图), 这部分内容先不放在这里。

weighted_unifrac.png-19.3kB

流程里有两点没有操作:

a. 最OTU表标准化
b. 去除线粒体和叶绿体来源的16S序列。

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

Last Update: 2017-11-01 3:20 PM

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

引言

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序列丰度。

三、DADA2实践

3.1 准备材料:

下载测试文件

wget http://www.mothur.org/w/images/d/d6/MiSeqSOPData.zip

下载训练集

Silva 128 [Silva taxonomic training data formatted for DADA2 (Silva version 128)][65]

wget https://zenodo.org/record/824551/files/silva_nr_v128_train_set.fa.gz
wget https://zenodo.org/record/824551/files/silva_species_assignment_v128.fa.gz

RDP trainset 16 [RDP taxonomic training data formatted for DADA2][76]

wget https://zenodo.org/record/801828/files/rdp_species_assignment_16.fa.gz
wget https://zenodo.org/record/801828/files/rdp_train_set_16.fa.gz

3.2 DADA2 安装

测试平台:R 3.3.0 wget https://zenodo.org/record/158955/files/rdp_species_assignment_14.fa.gz
wget https://zenodo.org/record/158955/files/rdp_train_set_14.fa.gz

UNITE
wget https://unite.ut.ee/sh_files/sh_general_release_10.10.2017.zip

DADA2 Pipeline

测试平台:R 3.3.0

DADA2 R Pakcage 安装

使用默认安装模式:

source("https://bioconductor.org/biocLite.R")
biocLite("dada2")

当前安装版本为: 1.2.2

使用devtools安装模式:

devtools::install_github("benjjneb/dada2")

安装最新版本: 1.5.8

3.3 测试流程

“`r
library(dada2);
packageVersion(“dada2”)

rdp_path <- "/project/task/dada2/RDP-16"
miseq_sop <- "/project/task/dada2/MiSeq_SOP"

fnFs <- sort(list.files(miseq_sop, pattern="_R1_001.fastq"))
fnRs <- sort(list.files(miseq_sop, pattern="_R2_001.fastq"))

sample.names <- sapply(strsplit(fnFs, "_"), `[`, 1)
fnFs <- file.path(miseq_sop, fnFs)
fnRs <- file.path(miseq_sop, fnRs)

plotQualityProfile(fnFs[1:2])

filt_path <- file.path(miseq_sop, "filtered")
filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))

filt <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
      maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
      compress=TRUE, multithread=TRUE)

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)

names(derepFs) <- sample.names
names(derepRs) <- sample.names

errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)

dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
seqtab  <- makeSequenceTable(mergers)
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)

getN <- function(x) sum(getUniques(x))
track <- cbind(filt, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled", "nonchim")
rownames(track) <- sample.names

taxa <- assignTaxonomy(seqtab.nochim,  file.path(rdp_path, "rdp_train_set_16.fa.gz"), multithread=TRUE)
seqtab_taxa <- cbind('#seq'=rownames(taxa), t(seqtab.nochim), taxa)
write.table(seqtab_taxa, "dada2-counts-taxon.txt", sep="\t", quote=F, row.names = F)

taxa.plus <- addSpecies(taxa, file.path(rdp_path, "rdp_species_assignment_16.fa.gz"))
seqtable.taxa.plus <- cbind('#seq'=rownames(taxa.plus), t(seqtab.nochim), taxa.plus)
write.table(seqtable.taxa.plus , "dada2-counts-taxon.species.txt", sep="\t", quote=F, row.names = F)

“`

3.4 流程解读

整个流程包含以下几个步骤:

序列过滤和切除
错误率计算
去噪音
合并双端序列
去除嵌合体
物种分类
物种水平鉴定

整个流程通过几个关键函数完成,filterAndTrim、derepFastq、learnErrors、dada、mergePairs、assignTaxonomy、addSpecies 等几个核心函数。

值得关注是addSpecies和assignSpecies, 逻辑是通过完全匹配的序列分类一致性(重点在这里, 两个问题: 1. 库应该足够大, 2. 序列库足够准确)去判断是否应该分配species信息;

DADA2 提供了silva_species_assignment_v128.fa.gzrdp_species_assignment_16.fa.gz用于识别可以分类到种水平信息, 该文件是通过对原始序列问题进行几个操作实现:

a. 去除分类不明确的序列,比如 Uncultured、 unclassified、Outgroup、Unidentified序列;
b. 去除分类不完整序列,要包含属和种水平信息;
c. 只保留属和种信息;

3.5 QIIME整合

QIIME 通过插件形式整合DADA2用于构建这中feature-table, 代替传统的OTU。

qiime dada2 命令可以执行特征表构建, 并使用MD5描述信息,支持单端(denoise-single)和双端模式(denoise-paired),

支持的命令行接口(参数基本可以和DADA2的R函数对应):

--i-demultiplexed-seqs: demultiplexe后的序列;
--p-trunc-len-f: 针对forward序列,截断序列长度大小,设置为0不执行截断操作;
--p-trunc-len-r: 针对reverse序列,截断序列长度大小,设置为0不执行截断操作;
--o-table: FeatureTable,用于下游各种分析。
--o-representative-sequences: 代表序列,簇心。
--p-hashed-feature-ids:  使用Hash值做序列标识符。
--p-no-hashed-feature-ids:不使用Hash值做序列标识符。
--p-n-reads-learn: 用于无监督学习的序列数。
--p-n-threads: 线程数;

实例:

qiime dada2 denoise-paired \
    --i-demultiplexed-seqs paired-end-sequences.qza \
    --p-trunc-len-f  0   \
    --p-trunc-len-r  0   \
    --o-table  table.qza \
    --o-representative-sequences \
    rep-seqs.qza

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

Last Upate: 2017/10/16 21:20

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

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

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

  1. The 97% threshold is much too low, especially for V4
  2. 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 both  : 支持反向互补序列;

实例:

usearch -fastx_uniques seq.fasta -fastaout uniques.fasta -sizeout -relabel Uniq
2. fastx_uniques_persample

fastx_uniques_persample 为 UsearchV10 版本新的函数,实现了标题需求。

支持的可选参数比较少, 见官方文档

-fastaout     : fasta文件格式输出;
-sizeout      : 输出序列条数信息;
-minuniquesize: 最小序列数/per sample 阈值,如果低于这个数值,序列不输出。

命令行模式用例:

usearch  -fastx_uniques_persample seq.fna  -fastaout  unique.fa  -sizeout -minuniquesize 1

输入文件格式要求: 序列必须使用 sample=xxx 对序列进行注释;

3. DADA2

DADA2 通过makeSequenceTable函数实现计算每个样本中Uniques 序列丰度谱。

derep1 <- derepFastq(system.file("extdata", "sam1F.fastq.gz", package="dada2"))
derep2 <- derepFastq(system.file("extdata", "sam2F.fastq.gz", package="dada2"))
dada1 <- dada(derep1, tperr1)
dada2 <- dada(derep2, tperr1)
makeSequenceTable(list(sample1=dada1, sample2=dada2))

生成类似OTU的数值矩阵,其实逻辑很简单,就是查数。

4. seqtk_seqtab

fastx_uniques_persample 输出格式只能是fasta文件格式,还需要后续处理,才能得到序列丰度表,所以就再次造了轮子: seqtk_seqtab, seqtk_seqtab 实现是因为DADA2的模式, 毕竟使用 R函数 没有直接使用命令行来的快。

命令接口:

$seqtk_seqtab
Usage: seqtk_seqtab <fastq>
version: 0.0.1

实例:

seqtk_seqtab seq.fna > seqtab.txt

速度很快,约为 fastx_uniques_persample 一半时间。

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

Last Update: 2017-10-27 4:17 PM

扩增子数据分析之标准化:抽样 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进行计算:

计算标准化需要将浮点数转换成整数,为了避免太多小数,比如0, USEARCH推荐大于2倍样本中典型样本大小, 比如大部分都是30000条, 这时可以选择 60000条。

计算公式:

new = (old * sample_size) / total 

实例:

usearch -otutab_norm otutable_counts.txt  -sample_size 50000  -output otutab_norm.txt
3. VST变换

phyloseq 使用DESeq / DESeq2 进行标准化,并计算丰度差异差异OTU, QIIME也提供了normalize_table.py 用于标准化原始OTU count数。
DESEQ2 提供了 varianceStabilizingTransformation函数对counts进行变化,即VST变换。

实例:

varianceStabilizingTransformation(dds, blind = TRUE, fitType = "parametric")
ps_dds <- estimateSizeFactors(ps_dds)
ps_dds <- estimateDispersions(ps_dds)
getVarianceStabilizedData(dds)

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

Last Update: 2017-10-31 3:36 PM

扩增子测序分析之多样性指数: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

即:

Bray-Curtis = 1 – 2*S / T

其中:

S = sum k = 1 .. N of min(n_ki, n_kj)
T = sum k = 1 .. N of (n_ki + n_kj)    
N = number of OTUs.
n_ki = count of OTU k in sample i.

如果不考虑数值因素, 可以即将OTU丰度大于1的数值设置成1, 为0的数值为0,这样表征为存在/缺失模式。

2.2 Jaccard distance

Jaccard distance 可通过Jaccard similarity coefficient获得,即:d= 1 - J

Jaccard = 1 – S / T

其中:

S = sum k = 1 .. N of min(n_ki, n_kj)
T = sum k = 1 .. N of max(n_ki, n_kj)
N = number of OTUs.
n_ki = count of OTU k in sample i.

如果不考虑数值因素, 可以即将OTU丰度大于1的数值设置成1, 为0的数值为0,这样表征为存在/缺失模式。

Jaccard distance

该参数被广泛应用到比较相似度度量,比如使用Kmer度量两个基因组有多相似来鉴定菌株。

2.3 UniFrac 首先使用OTU代表序列构建系统树,比如:

OTU-Tree

节点描述OTU,每个节点(包括叶子节点和中间节点)之间定义分支,图上标注的长度描述节点之间距离。

Quantitative and Qualitative β Diversity Measures Lead to Different Insights into Factors That Structure Microbial Communities 描述了如何计算 Unweighted UniFrac和weighted UniFrac距离,如下图所示:

UniFrac

在计算UniFrac距离时一般会先修剪系统进化树,去除不包含代比较样本的节点和分支。

Unweighted UniFrac距离:

Unweighted UniFrac定义为特异的分支长度/观测到的分支长度。

UniFrac.png

上图(a)中灰色的线为共享的分支。黑色分支为特异的分支。

比如下面几种情况:

STAMPS2015_Knight

图片来自 STAMPS2015_Knight 幻灯片。

Weighted UniFrac距离:

Weighted UniFrac计算当时比较复杂,需要使用样本相对丰度对分支进行加权: 首先需要原始weighted UniFrac value (u)值,通过下面公式:

u-value

其中: N : 为树中的所有分支; b_i : 分支i的长度; A_i和B_i: 样本A和B中来源分支i序列数; A_T和B_T: 样本A和B中的序列数;

3. 计算距离矩阵

3.1 USEARCH div_beta

USEARCH div_beta 提供了比较简单的接口计算距离矩阵。 支持的命令行参数有:

-filename_prefix  :指定文件输出路径;
‑mx_suffix        :矩阵前缀;
‑sorted_mx_suffix :排序的矩阵前缀(相似样本排在一起);
-tree_suffix      :生成的距离树前缀(等级聚类树);
-metrics          :指定多样性指数类型,支持bray_curtis,bray_curtis_binary,
                   euclidean,jaccard,jaccard_binary,manhatten,unifrac,unifrac_binary
-tree             :计算Unifrac需要提供系统进化树信息(Newick格式)

如果添加任何选项,则计算所有支持的距离矩阵。

usearch -beta_div otutable_counts.txt -metrics bray_curtis,jaccard -filename_prefix /biostack/task/otutable-utils/otus/beta_div/

也可以通过

usearch -beta_div otutable_counts.txt -metrics unifrac -tree  rep_set.nwk   -filename_prefix /biostack/task/otutable-utils/otus/beta_div/

3.2 QIIME2 div_beta

qiime2 数据分析流程通过 qiime diversity接口提供了分析beta多样的各种命令:

支持的beta多样性指数:

--i-table : FeatureTable
--p-n-jobs:并行参数;
--p-metric: correlation|yule|sqeuclidean|sokalsneath|seuclidean|braycurtis
            |canberra|hamming|matching|dice|wminkowski|cosine|sokalmichener
            |russellrao|chebyshev|rogerstanimoto|cityblock|jaccard
            |kulsinski|mahalanobis|euclidean
--o-distance-matrix: 输出距离矩阵;
--output-dir: 输出目录(如不指定--o-distance-matrix);

实例:

qiime diversity  beta        \
      --p-n-jobs 2           \
      --i-table table.qza    \
      --p-metric braycurtis  \
      --o-distance-matrix  distance-matrix.qza

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

Last Update: 2017-09-25 11:56 AM

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

1.稀释曲线介绍

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

Rarefaction Curve

2.usearch alpha_div_rare

Usearch 提供了alpha_div_rarecmd_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: 最大稀释深度,必须大于min_depth.
--p-min-depth:最小稀释深度, 默认1;
--p-steps: 稀释步长,默认10;
--p-iterations: 迭代次数 默认10;
--i-phylogeny: 如计算PDTree, 需要提供系统进化树;
--p-n-jobs:并行参数;
--p-metric: enspie|michaelis_menten_fit|faith_pd|lladser_pe|fisher_alpha
            |goods_coverage|doubles|simpson|margalef|observed_otus
            |shannon|pielou_e|chao1|brillouin_d|menhinick|simpson_e
            |robbins|dominance|heip_e|singles|mcintosh_d|ace
            |mcintosh_e|gini_index|berger_parker_d
--o-visualization: 输出稀释曲线;
--output-dir: 输出目录(如不指定--o-visualization);

实例:

qiime  diversity  alpha-rarefaction  \
     --i-table   table.qza      \
      --p-max-depth 1000        \
      --o-visualization  alpha-rarefaction.qza \
      --p-metrics  simpson      \
      --p-metrics observed_otus \
      --p-metrics shannon

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

Last Update: 2017-09-25 11:56 AM