Archives

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

Comments are closed.