qiime 2分析菌群数据代码

1. installation

有Linux服务器的伙伴推荐使用Conda安装，想在windows笔记本上体验的朋友可使用Virtualbox虚拟机安装并学习

#install conda
conda update conda  #升级conda程序
conda install wget  #安装下载工具

#install QIIME2 with conda
wget https://data.qiime2.org/distro/core/qiime2-2018.2-py35-linux-conda.yml  #下载软件安装列表
conda env create -n qiime2-2018.2 --file qiime2-2018.2-py35-linux-conda.yml  #创建虚拟环境并安装qiime2，防止影响其他已安装软件
# OPTIONAL CLEANUP
rm qiime2-2018.2-py35-linux-conda.yml  #删除软件列表

#activate qiime2
source activate qiime2-2017.12  #激活工作环境

#deactivate qiime2
souce deactivate   #关闭工作环境

2.Obtaining and importing data
    

#singleSingleEndFastqManifestPhred33
qiime tools import \
  --type 'SampleData[SequencesWithQuality]' \
  --input-path se-33-manifest \
  --output-path single-end-demux.qza \
  --source-format SingleEndFastqManifestPhred33

    #PairedEndFastqManifestPhred33
qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path pe-33-manifest \
  --output-path paired-end-demux.qza \
  --source-format PairedEndFastqManifestPhred33

qiime demux summarize \
  --i-data single-end-demux.qza \
  --o-visualization single-end-demux.qzv

3. Sequence quality control and feature table construction

qiime dada2 denoise-single \
  --i-demultiplexed-seqs single-end-demux.qza \
  --p-trim-left 0 \
  --p-trunc-len 165 \
  --o-representative-sequences rep-seqs-dada2.qza \
  --o-table table-dada2.qza \
  --o-denoising-stats stats-dada2.qza

mv rep-seqs-dada2.qza rep-seqs.qza
mv table-dada2.qza table.qza

4. FeatureTable and FeatureData summaries

qiime feature-table summarize \
  --i-table table.qza \
  --o-visualization table.qzv \
  --m-sample-metadata-file sample-metadata.tsv

#代表序列统计
qiime feature-table tabulate-seqs \
  --i-data rep-seqs.qza \
  --o-visualization rep-seqs.qzv

#qiime tools view rep-seqs.qzv # 展打开网页版统计结果

5. Generate a tree for phylogenetic diversity analyses

#多序列比对
qiime alignment mafft \
  --i-sequences rep-seqs.qza \
  --o-alignment aligned-rep-seqs.qza

#移除高变区
qiime alignment mask \
  --i-alignment aligned-rep-seqs.qza \
  --o-masked-alignment masked-aligned-rep-seqs.qza

#建树
qiime phylogeny fasttree \
  --i-alignment masked-aligned-rep-seqs.qza \
  --o-tree unrooted-tree.qza

#无根树转换为有根树
qiime phylogeny midpoint-root \
  --i-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza

6.Alpha and beta diversity analysis

#计算多样性
qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table table.qza \
  --p-sampling-depth 12052 \
  --m-metadata-file sample-metadata.tsv \
  --output-dir core-metrics-results

# 输出结果包括多种多样性结果，文件列表和解释如下：

# beta多样性bray_curtis距离矩阵 bray_curtis_distance_matrix.qza 

# alpha多样性evenness(均匀度，考虑物种和丰度)指数 evenness_vector.qza

# alpha多样性faith_pd(考虑物种间进化关系)指数 faith_pd_vector.qza

# beta多样性jaccard距离矩阵 jaccard_distance_matrix.qza

# alpha多样性observed_otus(OTU数量)指数 observed_otus_vector.qza

# alpha多样性香农熵(考虑物种和丰度)指数 shannon_vector.qza

# beta多样性unweighted_unifrac距离矩阵，不考虑丰度 unweighted_unifrac_distance_matrix.qza

# beta多样性unweighted_unifrac距离矩阵，考虑丰度 weighted_unifrac_distance_matrix.qza

# 统计faith_pd算法Alpha多样性组间差异是否显著，输入多样性值、实验设计，输出统计结果
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results/faith_pd_vector.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization core-metrics-results/faith-pd-group-significance.qzv

#统计evenness组间差异是否显著
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results/evenness_vector.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization core-metrics-results/evenness-group-significance.qzv

#网页展示结果，只要是qzv的文件，均可使用qiime tools view 查看或在线http://view.qiime2.org/查看
#qiime tools view evenness-group-significance.qzv


#按BodySite分组，统计unweighted_unifrace距离的组间是否有显著差异
qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-category BodySite \
  --o-visualization core-metrics-results/unweighted-unifrac-body-site-significance.qzv \
  --p-pairwise

# 可视化三维展示unweighted-unifrac的主坐标轴分析
qiime emperor plot \
  --i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \
  --m-metadata-file sample-metadata.tsv \
  --p-custom-axes DaysSinceExperimentStart \
  --o-visualization core-metrics-results/unweighted-unifrac-emperor-DaysSinceExperimentStart.qzv

#Alpha rarefaction plotting
qiime diversity alpha-rarefaction \
  --i-table table.qza \
  --i-phylogeny rooted-tree.qza \
  --p-max-depth 12052 \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization alpha-rarefaction.qzv

7.Training feature classifiers with q2-feature-classifier

wget \
  -O "85_otus.fasta" \
  "https://data.qiime2.org/2018.11/tutorials/training-feature-classifiers/85_otus.fasta"

wget \
  -O "85_otu_taxonomy.txt" \
  "https://data.qiime2.org/2018.11/tutorials/training-feature-classifiers/85_otu_taxonomy.txt"

qiime tools import \
  --type 'FeatureData[Sequence]' \
  --input-path 85_otus.fasta \
  --output-path 85_otus.qza

qiime tools import \
  --type 'FeatureData[Taxonomy]' \
  --source-format HeaderlessTSVTaxonomyFormat \
  --input-path 85_otu_taxonomy.txt \
  --output-path ref-taxonomy.qza

qiime feature-classifier extract-reads \
  --i-sequences 85_otus.qza \
  --p-f-primer CCAGACTCCTACGGGAGGCAG \
  --p-r-primer CGTATTACCGCGGCTGCTG \
  --o-reads ref-seqs.qza

qiime feature-classifier fit-classifier-naive-bayes \
  --i-reference-reads ref-seqs.qza \
  --i-reference-taxonomy ref-taxonomy.qza \
  --o-classifier classifier.qza

qiime feature-classifier classify-sklearn \
  --i-classifier classifier.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza

qiime metadata tabulate \
  --m-input-file taxonomy.qza \
  --o-visualization taxonomy.qzv

8. Taxonomic analysis

#物种分类柱状图
qiime taxa barplot \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization taxa-bar-plots.qzv

#Differential abundance testing with ANCOM
    #filter to gut samples
qiime feature-table filter-samples \
  --i-table table.qza \
  --m-metadata-file sample-metadata.tsv \
  --p-where "BodySite='gut'" \
  --o-filtered-table gut-table.qza

    #imputation method
qiime composition add-pseudocount \
  --i-table gut-table.qza \
  --o-composition-table comp-gut-table.qza

   9.Subject

qiime composition ancom \
  --i-table comp-gut-table.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-column Subject \
  --o-visualization ancom-Subject.qzv

    #specific taxonomic level
qiime taxa collapse \
  --i-table gut-table.qza \
  --i-taxonomy taxonomy.qza \
  --p-level 6 \
  --o-collapsed-table gut-table-l6.qza

qiime composition add-pseudocount \
  --i-table gut-table-l6.qza \
  --o-composition-table comp-gut-table-l6.qza

qiime composition ancom \
  --i-table comp-gut-table-l6.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-column Subject \
  --o-visualization l6-ancom-Subject.qzv

#shell
##check shell
echo $0
##Bach shell
source tab-qiime
10.Importing phylogenetic trees (unrooted)

qiime tools import \
  --input-path rep_set.tre \
  --output-path unrooted-tree.qza \
  --type 'Phylogeny[Unrooted]'

qiime phylogeny midpoint-root \
  --i-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza

11.Importing BIOM v1.0.0
 

qiime tools import \
  --input-path otu_table_mc2_w_tax_no_pynast_failures.biom \
  --type 'FeatureTable[Frequency]' \
  --source-format BIOMV100Format \
  --output-path feature-table-1.qza

12.Feature Table summarize
 

qiime feature-table summarize \
  --i-table feature-table-1.qza \
  --o-visualization table.qzv \
  --m-sample-metadata-file sample-metadata.tsv

13.Alpha and beta diversity analysis

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table feature-table-1.qza \
  --p-sampling-depth 51162 \
  --m-metadata-file sample-metadata.tsv \
  --output-dir core-metrics-results

##alpha-group-significance
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results/faith_pd_vector.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization core-metrics-results/faith-pd-group-significance.qzv

qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results/evenness_vector.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization core-metrics-results/evenness-group-significance.qzv

14.beta-group-significance

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-column BodySite \
  --o-visualization core-metrics-results/unweighted-unifrac-body-site-significance.qzv \
  --p-pairwise

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-column Subject \
  --o-visualization core-metrics-results/unweighted-unifrac-subject-group-significance.qzv \
  --p-pairwise

##PCoA plot
qiime emperor plot \
  --i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \
  --m-metadata-file sample-metadata.tsv \
  --p-custom-axes DaysSinceExperimentStart \
  --o-visualization core-metrics-results/unweighted-unifrac-emperor-DaysSinceExperimentStart.qzv

qiime emperor plot \
  --i-pcoa core-metrics-results/bray_curtis_pcoa_results.qza \
  --m-metadata-file sample-metadata.tsv \
  --p-custom-axes DaysSinceExperimentStart \
  --o-visualization core-metrics-results/bray-curtis-emperor-DaysSinceExperimentStart.qzv

#Alpha rarefaction plotting
qiime diversity alpha-rarefaction \
  --i-table feature-table-1.qza \
  --i-phylogeny rooted-tree.qza \
  --p-max-depth 51162 \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization alpha-rarefaction.qzv

#Taxonomic analysis

#Differential abundance testing with ANCOM