16S多样性分析docker镜像发布及使用
16S多样性分析docker镜像发布及使用
众所周知微生物多样性分析用到的软件众多，对于初学者而言安装软件是一件非常麻烦的事情，这里我们组学大讲堂发布了专门用于微生物扩增子分析docker镜像，里面内置了几乎所有要用到的关于多样性分析的相关软件，方便大家使用，只需安装docker并下载该镜像即可，避免繁琐的软件安装问题；
扩增子分析流程介绍
目前主流的扩增子分析流程包括，qiime1，qiime2 ，mothur ，usearch等。qiime以其方便好用的流程，被广大研究人员官方使用。逐渐的qiime2会替代qiime1，但是由于数据分析周期较长很多人还在使用qiime1并且也习惯使用qiime1的分析流程。因此我们提供这个ampliseq1镜像仍然安装qiime1，但也会安装：mothur ，usearch 软件方便大家分析数据，qiime2会安装在另一个单独的镜像中，后续发布；
镜像下载使用方法
1）首先系统中docker安装与使用方法见课程：
2）镜像搜索、下载、启动使用方法：
> docker search omicsclass #搜索组学大讲堂提供的所有镜像
NAME DESCRIPTION STARS OFFICIAL AUTOMATED
omicsclass/gene-family gene-family analysis docker image 5
omicsclass/rnaseq RNA-seq analysis docker image build by omics… 3
omicsclass/gsds-v2 GSDS 2.0 – Gene Structure Display Server 1
omicsclass/reseq whole genome resequence analysis 1
omicsclass/biocontainer-base Biocontainers base Image centos7 1
omicsclass/blast-plus blast+ v2.9.0 0
omicsclass/isoseq3 isoseq3 v3.3.0 build by omicsclass 0
omicsclass/bwa BWA v0.7.17 build by omicsclass 0
omicsclass/blastall legacy blastall v2.2.26 0
omicsclass/sratoolkit SRAtoolkit v2.10.3 and aspera v3.9.9.177872 0
omicsclass/ampliseq-q2 Amplicon sequencing qiime2 v2020.2 image 0
omicsclass/ampliseq-q1 Amplicon sequencing qiime1 v1.9.1 image 0
omicsclass/samtools samtools v1.10 build by omicsclass 0
omicsclass/bsaseq NGS Bulk Segregant Analysis image 0
omicsclass/gwas gwas analysis images 0
> docker pull omicsclass/ampliseq-q1 #下载扩增子镜像
>docker run --rm -it -m 4G --cpus 1 -v D:\qiime1-16s:/work omicsclass/ampliseq-q1:latest #启动并进入镜像
安装软件及版本介绍
OTU分析相关软件
qiime1 v1.9.1
mothur v.1.25.0
usearch 10.0.240
usearch61 6.1.544
vsearch v2.15.0
序列处理相关分析包
flash v2.15.0 序列合并
bioperl
biopython
fastqc
multiqc
fastp
qiime分析相关的包
blast-2.2.22
blat-36
cdhit-3.1
muscle-3.8.31
rdpclassifier-2.2
uclust
功能注释相关分析软件
picrust-1.1.4
bugbase
数据分析与可视化相关的R包：
ggplot2
Ternary
DESeq2
edgeR
ggtree
vegan
pheatmap
高级分析相关包
randomForest 机器学习
scikit-learn 机器学习
circos 圈图绘制
Krona 物种丰度圈图
LefSe 差异比较分析
注意：以上软件包的列表只是部分举例，实际安装的包还有更多。
扩增子镜像使用案例
1.数据准备
fastmap.txt文件, 测序数据文件放在data文件夹中，物种注释数据库文件Greengene silva unite等放在database目录，需要这些测试数据的同学可以关注组学大讲堂公众号，并回复16s，即可得到；准备完成之后，目录结构如下：
[root@aefe86d682b1 13:58:49 /work/amplicon_demo]# tree
data
|-- ERR3975186_1.fastq.gz
|-- ERR3975186_2.fastq.gz
|-- ERR3975187_1.fastq.gz
`-- ERR3975187_2.fastq.gz
fastmap.txt
设置一些环境变量方便后续调用：
dbdir=/work/database
workdir=/work/amplicon_demo
datadir=$workdir/data
fastmap=$workdir/fastmap.txt
mkdir /work/tmp
export TMPDIR=/work/tmp #防止临时目录存储 不够
######################database
#各大数据库地址：
silva_16S_97_seq=$dbdir/SILVA_132_QIIME_release/rep_set/rep_set_16S_only/97/silva_132_97_16S.fna
silva_16S_97_tax=$dbdir/SILVA_132_QIIME_release/taxonomy/16S_only/97/taxonomy_7_levels.txt
greengene_16S_97_seq=$dbdir/gg_13_8_otus/rep_set/97_otus.fasta
greengene_16S_97_tax=$dbdir/gg_13_8_otus/taxonomy/97_otu_taxonomy.txt
silva_18S_97_seq=$dbdir/SILVA_132_QIIME_release/rep_set/rep_set_18S_only/97/silva_132_97_18S.fna
silva_18S_97_tax=$dbdir/SILVA_132_QIIME_release/taxonomy/18S_only/97/taxonomy_7_levels.txt
unite_ITS_97_seq=$dbdir/unite_ITS8.2/sh_refs_qiime_ver8_97_04.02.2020.fasta
unite_ITS_97_tax=$dbdir/unite_ITS8.2/sh_taxonomy_qiime_ver8_97_04.02.2020.txt
2.合并双端reads，flash
cd $workdir #回到工作目录
mkdir 1.merge_pe
for i in `cat $fastmap |grep -v '#'|cut -f 1` ;do
echo "RUN CMD: flash $datadir/${i}_1.fastq.gz $datadir/${i}_2.fastq.gz \
-m 10 -x 0.2 -p 33 -t 1 \
-o $i -d 1.merge_pe"
flash $datadir/${i}_1.fastq.gz $datadir/${i}_2.fastq.gz \
-m 10 -x 0.2 -p 33 -t 1 \
-o $i -d 1.merge_pe
done
3.对原始数据进行fastqc质控
cd $workdir #回到工作目录
mkdir 2.fastqc
#fastqc查看数据质量分布等
fastqc -t 2 $workdir/1.merge_pe/*extendedFrags.fastq -o $workdir/2.fastqc
#质控结果汇总
cd $workdir/2.fastqc
multiqc .
4.数据质控：对原始序列进行去接头，引物，删除低质量的reads等等
cd $workdir #回到工作目录
mkdir 3.data_qc
cd 3.data_qc
#利用fastp工具去除adapter
#--qualified_quality_phred the quality value that a base is qualified.
# Default 15 means phred quality >=Q15 is qualified. (int [=15])
#--unqualified_percent_limit how many percents of bases are allowed to be unqualified
#--n_base_limit if one read's number of N base is >n_base_limit,
# then this read/pair is discarded
#--detect_adapter_for_pe 接头序列未知 可设置软件自动识别常见接头
#
for i in `cat $fastmap |grep -v '#'|cut -f 1`; do
echo "RUN CMD: fastp --thread 1 --qualified_quality_phred 10 \
--unqualified_percent_limit 50 \
--n_base_limit 10 \
--length_required 300 \
--trim_front1 29 \
--trim_tail1 18 \
-i $workdir/1.merge_pe/${i}.extendedFrags.fastq \
-o ${i}.clean_tags.fq.gz \
--adapter_fasta $workdir/data/illumina_multiplex.fa -h ${i}.html -j ${i}.json"
fastp --thread 1 --qualified_quality_phred 10 \
--unqualified_percent_limit 50 \
--n_base_limit 10 \
--length_required 300 \
--trim_front1 29 \
--trim_tail1 18 \
-i $workdir/1.merge_pe/${i}.extendedFrags.fastq \
-o ${i}.clean_tags.fq.gz \
--detect_adapter_for_pe -h ${i}.html -j ${i}.json
done
去除嵌合体
cd $workdir #回到工作目录mkdir 4.remove_chimerascd 4.remove_chimeras
#去除嵌合体
for i in `cat $fastmap |grep -v '#'|cut -f 1`; do #相同重复序列合并 vsearch --derep_fulllength $workdir/3.data_qc/${i}.clean_tags.fq.gz \ --sizeout --output ${i}.derep.fa #去嵌合体 vsearch --uchime3_deno ${i}.derep.fa \ --sizein --sizeout \ --nonchimeras ${i}.denovo.nonchimeras.rep.fa #相同序列还原为多个 vsearch --rereplicate ${i}.denovo.nonchimeras.rep.fa --output ${i}.denovo.nonchimeras.fadone
#根据参考序列去除嵌合体for i in `cat $fastmap |grep -v '#'|cut -f 1`; do vsearch --uchime_ref ${i}.denovo.nonchimeras.fa \ --db $dbdir/rdp_gold.fa \ --sizein --sizeout --fasta_width 0 \ --nonchimeras ${i}.ref.nonchimeras.fadone
qiime1分析pick otu 聚类方法
cd $workdir #回到工作目录
mkdir 5.pick_otu_qiime
cd 5.pick_otu_qiime
#合并fasta文件，并加序列号
for i in `cat $fastmap |grep -v '#'|cut -f 1`; do
rename_fa_id.pl -f $workdir/4.remove_chimeras/$i.ref.nonchimeras.fa \
-n $i -out $i.fa
done
#合并fa文件到qiime.fasta 之后删除所有单个样本的fa文件
cat *fa >qiime.fasta
rm -f *fa
###方法1：pick_de_novo_otus.py
###输出qiime pick otu 参数，更多：http://qiime.org/scripts/pick_otus.html
echo pick_otus:denovo_otu_id_prefix OTU >> otu_params_de_novo.txt
echo pick_otus:similarity 0.97 >> otu_params_de_novo.txt
echo pick_otus:otu_picking_method uclust >> otu_params_de_novo.txt #sortmerna, mothur, trie, uclust_ref, usearch, usearch_ref, blast, usearch61, usearch61_ref,sumaclust, swarm, prefix_suffix, cdhit, uclust.
echo assign_taxonomy:reference_seqs_fp $silva_16S_97_seq >> otu_params_de_novo.txt
echo assign_taxonomy:id_to_taxonomy_fp $silva_16S_97_tax >> otu_params_de_novo.txt
echo assign_taxonomy:similarity 0.8 >>otu_params_de_novo.txt
echo assign_taxonomy:assignment_method uclust >>otu_params_de_novo.txt # rdp, blast,rtax, mothur, uclust, sortmerna如果是ITS/18S数据，建议数据库更改为UNITE，方法改为blast。详细使用说明，请读官方文档：http://qiime.org/scripts/assign_taxonomy.html
pick_de_novo_otus.py -i qiime.fasta -f -o pick_de_novo_otus -p otu_params_de_novo.txt
alpha_diversity 分析
cd $workdir #回到工作目录
mkdir 8.alpha_diversity
cd 8.alpha_diversity
#alpha多样性指数展示
biom summarize-table -i $workdir/5.pick_otu_qiime/pick_de_novo_otus/otu_table_clean_rare.biom
echo alpha_diversity:metrics observed_species,PD_whole_tree,shannon,chao1,simpson,goods_coverage > alpha_params.txt
alpha_rarefaction.py -f -i $workdir/5.pick_otu_qiime/pick_de_novo_otus/otu_table_clean.biom -m $fastmap -o ./ -p alpha_params.txt -t $workdir/5.pick_otu_qiime/pick_de_novo_otus/rep_set.tre --retain_intermediate_files --min_rare_depth 40 --max_rare_depth 2032 --num_steps 10
#多样性指数差异比较qiime 自带检验与绘图
compare_alpha_diversity.py \
-i alpha_div_collated/chao1.txt \
-o alpha_chao1_stats \
-m $fastmap \
-t nonparametric \
-c city
compare_alpha_diversity.py \
-i alpha_div_collated/chao1.txt \
-o alpha_chao1_stats \
-m $fastmap \
-t parametric \
-c city
beta_diversity 分析
cd $workdir #回到工作目录
mkdir 9.beta_diversity
cd 9.beta_diversity
echo beta_diversity:metrics binary_jaccard,bray_curtis,unweighted_unifrac,weighted_unifrac,binary_euclidean > beta_params.txt
#-e 设置抽平数
beta_diversity_through_plots.py -f -i $workdir/5.pick_otu_qiime/pick_de_novo_otus/otu_table_clean.biom -m $fastmap -o ./ -t $workdir/5.pick_otu_qiime/pick_de_novo_otus/rep_set.tre -e 2844 -p beta_params.txt
#beta多样性adonis检验
compare_categories.py --method adonis -i unweighted_unifrac_dm.txt -m $fastmap -c Treatment -o adonis_out -n 999
部分分析结果展示：
课程推荐：微生物扩增子分析课程实操 微生物16S/ITS/18S分析原理及结果解读
生物信息入门到精通必修基础课：linux系统使用、docker搭建生物信息分析环境、实验室linux生信分析平台搭建、linux命令处理生物大数据、perl入门到精通、perl语言高级、R语言画图、R语言快速入门与提高、python语言入门到精通
- 发表于 2020-09-02 13:58
- 阅读 ( 6180 )
- 分类：宏基因组
