Amplicon sequencing analysis pipeline through qiime2 platform

qiime2是擴增子數(shù)據(jù)分析的最佳平臺之一众眨，其提供了大量從原始data到統(tǒng)計分析的插件缀程，尤其是它的可重復(fù)分析且可擴展插件的理念使得其成為擴增子分析首選的平臺。

Platform

qiime2是擴增子數(shù)據(jù)分析的最佳平臺之一盲赊，其提供了大量從原始data到統(tǒng)計分析的插件，尤其是它的可重復(fù)分析且可擴展插件的理念使得其成為擴增子分析首選的平臺敷扫。對于如何安裝該平臺哀蘑，個人建議使用conda安裝，并且官網(wǎng)也提供了conda安裝的yaml文件葵第，但為了快速安裝绘迁，我們需要把conda鏡像重新設(shè)置修改一下。

step1: download the yaml file

wget https://data.qiime2.org/distro/core/qiime2-2020.8-py36-linux-conda.yml

step2: update the conda version

conda update conda -y

step3: modify the .condarc and yaml file

# reset the conda channels 
channels:
  - conda-forge
  - bioconda
  - biobakery
  - qiime2
  - ohmeta
  - defaults
show_channel_urls: true
channel_alias: https://mirrors.bfsu.edu.cn/anaconda
default_channels:
  - https://mirrors.bfsu.edu.cn/anaconda/pkgs/main
  - https://mirrors.bfsu.edu.cn/anaconda/pkgs/free
  - https://mirrors.bfsu.edu.cn/anaconda/pkgs/r
custom_channels:
  conda-forge: https://mirrors.bfsu.edu.cn/anaconda/cloud
  bioconda: https://mirrors.bfsu.edu.cn/anaconda/cloud
  biobakery: https://mirrors.bfsu.edu.cn/anaconda/cloud
  qiime2: https://mirrors.bfsu.edu.cn/anaconda/cloud
  ohmeta: https://mirrors.bfsu.edu.cn/anaconda/cloud
  knights-lab: https://conda.anaconda.org
  alienzj: https://conda.anaconda.org
  
 # change the yaml channel
 - qiime2/label/r2020.8
 + qiime2

step4: create the qiime2-2020.8 conda environment

conda env create -n qiime2-2020.8 --file qiime2-2020.8-py36-linux-conda_update.yml -y

step5: activate environment and test installation

# actiavate 
conda activate qiime2-2020.8 
# test 
qiime --help
# deactivate 
conda deactivate 

# Note: if you have trouble with no activate command, you need use *eval "$(conda shell.bash hook)"*

Analysis pipeline

Check the fastq phred score

通過fastqc軟件處理fastq獲取原始數(shù)據(jù)的reads質(zhì)量分布等情況卒密，為后續(xù)設(shè)置堿基質(zhì)量閾值做好準(zhǔn)備缀台。

fastqc --noextract -f fastq input.fq.gz  -o result/outdir

Get the positions of primer in the fastq

DADA2算法需要提供forward和reverse引物序列在reads上的位置信息，我們可以通過使用usearch和vsearch軟件獲取引物的位置信息栅受，并將其作為參數(shù)配置給qiime2 dada2插件将硝。獲取primer_hits.txt文件的開始和結(jié)束10個reads的primer信息恭朗，找到引物的位置信息。

# step1: merge fq
echo "step1: mkdir result and merger all PE fq files"
mkdir result
./usearch -fastq_mergepairs rawdata/*_R1.fq -fastqout result/all_samples_merged.fq -relabel @

# step2: sampling sequence
echo "step2: sampling sequence for primer check"
./usearch -fastx_subsample result/all_samples_merged.fq -sample_size 20000 -fastqout result/all_sub_for_primer_check.fq

# step3: search primer position
echo "step3: search primer position of the sequence"
./usearch -search_oligodb result/all_sub_for_primer_check.fq -db primers.fasta -strand both -userout result/primer_hits.txt -userfields query+qlo+qhi+qstrand

# step4: obtain the position information of primer
head primer_hits.txt  # head 23  tail 19
#HTN1.182        452     471     -
#HTN1.182        7       23      +
#HTN1.183        447     466     -
#HTN1.183        7       23      +
#HTN1.233        428     447     -
#HTN1.233        7       23      +
#HTN1.570        430     449     -
#HTN1.570        7       23      +
#HTN1.219        448     467     -
#HTN1.219        7       23      +
tail primer_hits.txt # head 23  tail 19
#Normal6.996134  429     448     -
#Normal6.996134  7       23      +
#Normal6.996331  430     449     -
#Normal6.996331  7       23      +
#Normal6.996257  448     467     -
#Normal6.996257  7       23      +
#Normal6.996360  429     448     -
#Normal6.996360  7       23      +
#Normal6.961965  430     449     -
#Normal6.961965  7       23      +

Convert fastq into qza

step1: 準(zhǔn)備符合import格式的文件 manifest.tsv(包含樣本ID以及forward和reverse fq路徑的文件)和sample-metadata.tsv(樣本的分組信息)

# generate manifest.tsv 
find /rawdata/ -name "*gz" | grep '1_' | sort | perl -e 'print"sampleid\tforward-absolute-filepath\treverse-absolute-filepath\n"; while(<>){chomp; $name=(split("\/", $_))[-1]; $name1=$name; $name2=$_; $name1=~s/_[1|2]_.fq.gz//g; $name2=~s/1_/2_/; print "$name1\t$_\t$name2\n";}'  > pe-33-manifest-trimmed.tsv

# sample-metadata.tsv details
#sampleid        Treatment
#HTN1    HTN
#HTN2    HTN
#Normal1 Normal

step2: import fastq into qza依疼，雙端和單端數(shù)據(jù)的import參數(shù)不同痰腮，并且不同的phred score使用的參數(shù)也不一樣

# PE mode 
qiime tools import \
      --type 'SampleData[PairedEndSequencesWithQuality]' \
      --input-path pe-33-manifest-trimmed.tsv \
      --output-path result/paired-end-demux.qza \
      --input-format PairedEndFastqManifestPhred33V2

# SE mode 
qiime tools import \
  --type "SampleData[SequencesWithQuality]" \
  --input-path single-33-manifest.tsv \
  --output-path result/single-end-demux.qza \
  --input-format SingleEndFastqManifestPhred33V2 

Sequence quality control

Step1: 為了更加準(zhǔn)確地過濾低質(zhì)量的堿基，可以再使用qiime2自帶summarize插件查看低質(zhì)量堿基的位置分布律罢，最后再結(jié)合第二步usearch和vsearch的primer位置信息設(shè)置適合過濾的參數(shù)膀值。

qiime demux summarize \
  --i-data result/paired-end-demux.qza \
  --o-visualization result/paired-end-demux-summary.qzv

Step2: DADA2算法相比常用的OTU算法，其計算的amplicon variant sequences（ASV）的feature會更好一些误辑，feature代替OTU是一種趨勢沧踏。在此之后，Usearch的開發(fā)者Robert C Edgar迅速開發(fā)了更好的unoise2算法巾钉，該算法已更新到unoise3翘狱，并放話unoise2比DADA2更準(zhǔn)確。

DADA2是R的一個軟件包砰苍，可以進(jìn)行過濾潦匈，去重，嵌合體過濾赚导，reads的拼接茬缩，可以修正擴增子的測序錯誤，確定更多的真實變異吼旧。擴增子測序本身就具有內(nèi)在的限制凰锡，但是聚類OTU的方式進(jìn)一步限制了它的發(fā)展。OTU不是物種圈暗，它們不應(yīng)該成為錯誤的一部分掂为，DADA2可以具有更高的分辨率

DADA(Divisive Amplicon Denoising Algorithm)含義為區(qū)分?jǐn)U增子降噪方程可以確定真實的變異在454測序擴增子數(shù)據(jù)輸出更少的假陽性。DADA2是DADA的擴展和增強可以應(yīng)用于Illumina測序數(shù)據(jù)

• 特點：DADA2最重要的優(yōu)勢是它用了更多的數(shù)據(jù)厂置。DADA2的錯誤模型包含了質(zhì)量信息菩掏，而其他的方法都在過濾低質(zhì)量之后把序列的質(zhì)量信息忽略。而且DADA2的錯誤模型也包括了定量的豐度昵济，而且該模型也計算了各種不同轉(zhuǎn)置的概率A->C智绸。而且DADA2以自身數(shù)據(jù)的錯誤模型為參數(shù)，不用依賴于其他參數(shù)分布模型访忿。
  DADA2算法：一種分列式算法
• 原理：
  • 1 首先將每個reads全部看作單獨的單元瞧栗，Sequence相同的reads被納入一個sequence，reads個數(shù)即成為該sequence的豐度（abundance）（其實就是去冗余的過程）
  • 2 計算每個sequence豐度的p-value海铆。當(dāng)最小的p-value低于設(shè)定的閾值時迹恐， 將產(chǎn)生一個新的partition。每一個sequence將會被歸入最可能生成該sequence的partition卧斟。
  • 3 依次類推殴边，完成分割歸并憎茂。

# DADA2 denosie
qiime dada2 denoise-paired \
        --i-demultiplexed-seqs result/paired-end-demux.qza \
        --p-trim-left-f 23 \
        --p-trim-left-r 19 \
        --p-trunc-len-f 0 \
        --p-trunc-len-r 0 \
        --p-n-threads 20 \
        --o-table result/table.qza \
        --o-representative-sequences result/rep-seqs.qza \
        --o-denoising-stats result/stats.qza

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

Taxonomic annotation

step1: 通過比對已知分類學(xué)組成的參考數(shù)據(jù)庫的序列，可以獲知feature table的代表序列的物種注釋情況锤岸。在qiime2通呈＃可以使用已經(jīng)搭建好的分類學(xué)分類器：silva132和Greengene 13_8等。

GreenGene數(shù)據(jù)庫比較明顯的問題就是屬種水平注釋低是偷，所以很多條目里拳氢，g和s下劃線后面都是空的，如果關(guān)注屬種水平的注釋蛋铆，則不建議使用該數(shù)據(jù)庫馋评。

結(jié)合相關(guān)專業(yè)人士的反饋意見，個人建議使用silva數(shù)據(jù)庫作為物種注釋的首選參考數(shù)據(jù)庫刺啦。

# downlaod silva classifier data 
wget https://data.qiime2.org/2020.8/common/silva-138-99-nb-classifier.qza 

# annotation
qiime feature-classifier classify-sklearn \
                --i-classifier database/silva-138-99-nb-classifier.qza  \
                --i-reads result/rep-seqs.qza \
                --o-classification result/taxonomy-dada2-sliva.qza \
                --p-n-jobs 20 \
                --verbose \
                --output-dir result/

step2: 可通過將某些代表序列與擴增子數(shù)據(jù)庫在blast軟件下再進(jìn)行物種注釋留特，該結(jié)果與qiime2提供的分類學(xué)分類器結(jié)果比較，從而可以評估分類學(xué)分類器的性能玛瘸，這適合在構(gòu)造新的分類學(xué)分類器時候使用磕秤。

# extract the validated sequences by qiime2 view network site
qiime feature-table tabulate-seqs \
   --i-data result/rep-seqs.qza \
   --o-visualization result/rep-seqs.qzv

Filter the unsuitable ASV

step1: remove low occurrence ASVs

根據(jù)table的結(jié)果設(shè)置過濾threshold，閾值有frequency和samples捧韵，即ASV在所有樣本的總reads和出現(xiàn)在樣本數(shù)目。計算平均采樣深度（對所有ASV的count加和并求平均值）汉操，設(shè)置采樣閾值后再乘以平均采樣深度即獲得frequency閾值再来，另外也可以設(shè)置ASV出現(xiàn)在多少樣本內(nèi)。

qiime feature-table filter-features \
   --i-table result/table.qza \
   --p-min-frequency 10 \
   --p-min-samples 1 \
   --o-filtered-table result/table_filter_low_freq.qza

step2: remove contamination and mitochondria, chloroplast sequence.

16s擴增子常見污染序列是來自于線粒體和葉綠體等16s序列磷瘤，另外也存在一些未注釋的序列芒篷，均需要去除。

qiime taxa filter-table \
   --i-table result/table_filter_low_freq.qza \
   --i-taxonomy result/taxonomy-dada2-sliva.qza \
   --p-exclude mitochondria,chloroplast \
   --o-filtered-table  result/table_filter_low_freq_contam.qza 

step3: drop the low depth samples:

經(jīng)過上述處理后采缚，某些樣本含有較少的ASV總量针炉，因此可以將其剔除。通常使用的threshold的范圍是1,000 - 4,000 reads扳抽。

# summarise all the ASV counts in each sample 
qiime feature-table summarize \
   --i-table result/table_filter_low_freq_contam.qza \
   --o-visualization result/table_filter_low_freq_contam_summary.qzv
# remove samples
qiime feature-table filter-samples \
   --i-table  result/table_filter_low_freq_contam.qza \
   --p-min-frequency 4000 \
   --o-filtered-table  result/final_table.qza 
# representative sequence
qiime feature-table filter-seqs \
   --i-data result/rep-seqs.qza \
   --i-table result/final_table.qza \
   --o-filtered-data result/final_rep_seqs.qza
# reannotate
qiime feature-classifier classify-sklearn \
   --i-classifier database/silva-138-99-nb-classifier.qza  \
   --i-reads result/final_rep_seqs.qza \
   --o-classification result/final_taxonomy_sliva.qza \
   --p-n-jobs 20 \
   --verbose \
   --output-dir result/
# core features
qiime feature-table core-features \
   --i-table result/final_table.qza \
   --p-min-fraction 0.6 \
   --p-max-fraction 1 \
   --p-steps 11 \
   --o-visualization result/final_table_cores.qzv \
   --output-dir result

Downstream analysis

Constructing phylogenetic tree and diversity analysis

step1: 系統(tǒng)發(fā)育樹能夠服務(wù)于后續(xù)多樣性分析

qiime phylogeny align-to-tree-mafft-fasttree \
      --i-sequences result/final_rep_seqs.qza \
      --o-alignment result/final_rep_seqs_aligned.qza \
      --o-masked-alignment result/final_rep_seqs_masked.qza \
      --p-n-threads 20 \
      --o-tree result/unrooted-tree.qza \
      --o-rooted-tree result/rooted-tree.qza

step2: rarefication curve:

稀疏曲線可以了解測序深度與ASV的關(guān)系

qiime diversity alpha-rarefaction \
     --i-table result/final_table.qza \
     --i-phylogeny result/rooted-tree.qza \
     --p-max-depth 60000 \
     --m-metadata-file sample-metadata.tsv \
     --o-visualization result/p-max-depth-60000-alpha-rarefaction.qzv

step3: diversity analysis

根據(jù)ASV的最小測序深度設(shè)置sampling參數(shù)

# all diversity index and distance 
qiime diversity core-metrics-phylogenetic \
     --i-phylogeny result/rooted-tree.qza \
     --i-table result/final_table.qza \
     --p-sampling-depth 60000 \
     --m-metadata-file sample-metadata.tsv \
     --output-dir result/sample-depth-60000-core-metrics-results    

step4: faith_pd diversity parameters

# example for faith_pd_vector of group analysis
qiime diversity alpha-group-significance \
     --i-alpha-diversity result/sample-depth-60000-core-metrics-results/faith_pd_vector.qza \
     --m-metadata-file sample-metadata.tsv \
     --o-visualization result/sample-depth-60000-core-metrics-results/faith-pd-group-significance.qzv
# example for alpha diversity of group analysis
qiime diversity alpha-group-significance \
     --i-alpha-diversity result/sample-depth-60000-core-metrics-results/shannon_vector.qza \
     --m-metadata-file sample-metadata.tsv \
     --o-visualization result/shannon_compare_groups.qzv 
# beta diversity 
qiime diversity beta-group-significance \
    --i-distance-matrix result/sample-depth-60000-core-metrics-results/unweighted_unifrac_distance_matrix.qza \
    --m-metadata-file sample-metadata.tsv \
    --m-metadata-column Treatment \
    --p-pairwise false \
    --p-permutations 999 \
    --o-visualization result/unweighted-unifrac-subject-significance.qzv 
# three dimensions to show beta diversity
qiime emperor plot \
    --i-pcoa result/sample-depth-60000-core-metrics-results/unweighted_unifrac_pcoa_results.qza \
    --m-metadata-file sample-metadata.tsv \
    --p-custom-axes Treatment \
    --o-visualization result/unweighted-unifrac-emperor-height.qzv

visualizing taxonomic composition

qiime taxa barplot \
     --i-table result/final_table.qza \
     --i-taxonomy result/final_taxonomy_sliva.qza \
     --m-metadata-file sample-metadata.tsv \
     --o-visualization result/final_taxa_barplots_sliva.qzv

Analysis of composition of microbiomes (ANCOM)

ANCOM（可以了解下sparse compositional correlation (SparCC) to analyze correlation networks among taxa）可用于比較微生物在組間差異的分析方法篡帕， 結(jié)果與LEfse類似。該方法基于成分對數(shù)比的方法贸呢，即先對count數(shù)據(jù)進(jìn)行對數(shù)轉(zhuǎn)換镰烧，再通過簡單的秩和檢驗（stats包內(nèi)的aov, friedman.test, lme等函數(shù)）進(jìn)行比較，最后計算統(tǒng)計量w楞陷。ANCOM的結(jié)果用W值來衡量組間差異顯著性怔鳖。W值越高代表該物種在組間的差異顯著性越高。ANCOM的R代碼(推薦９潭辍＝嶂础度陆！)。

# add pseudocount for log transform
qiime composition add-pseudocount \
   --i-table result/final_table.qza \
   --p-pseudocount 1 \
   --o-composition-table result/final_table_pseudocount.qza
# ANCOM 
qiime composition ancom \
   --i-table result/final_table_pseudocount.qza \
   --m-metadata-file sample-metadata.tsv \
   --m-metadata-column Treatment \
   --output-dir result/ancom_output

export qza into other format type data

qza數(shù)據(jù)文件

QIIME2為了使分析流程標(biāo)準(zhǔn)化献幔，分析過程可重復(fù)懂傀，制定了統(tǒng)一的分析過程文件格式.qza；qza文件類似于一個封閉的系統(tǒng)斜姥，里面包括原始數(shù)據(jù)鸿竖、分析的過程和結(jié)果；這樣保證了文件格式的標(biāo)準(zhǔn)铸敏，同時可以追溯每一步的分析缚忧，以及圖表繪制參數(shù)。這一方案為實現(xiàn)將來可重復(fù)的分析提供了基礎(chǔ)杈笔。

# representative sequences
qiime tools export \
   --input-path result/final_rep_seqs.qza \
   --output-path final_result
# features table
qiime tools export \
   --input-path result/final_table.qza \
   --output-path final_result
biom normalize-table \
    -i final_result/feature-table.biom \
    -r \
    -o final_result/feature-table-norm.biom
biom convert \
    -i final_result/feature-table-norm.biom \
    -o final_result/feature-table-norm.tsv \
    --to-tsv \
    --header-key taxonomy

LEfse

LEfse是LDA Effect Size分析闪水，其本質(zhì)是一類判別分析。其結(jié)果一般配合進(jìn)化分支圖使用蒙具，也即是展示差異物種在進(jìn)化上的關(guān)系球榆。推薦使用yintools的LEfse的R腳本。remotes::install_github("ying14/yingtools2")

原理：首先使用non-parametric factorial Kruskal-Wallis (KW) sum-rank test（非參數(shù)因子克魯斯卡爾—沃利斯和秩驗檢）檢測具有顯著豐度差異特征禁筏，并找到與豐度有顯著性差異的類群持钉。最后，LEfSe采用線性判別分析（LDA）來估算每個組分（物種）豐度對差異效果影響的大小篱昔。

進(jìn)化分支圖：由內(nèi)至外輻射的圓圈代表了由門至屬（或種）的分類級別每强。在不同分類級別上的每一個小圓圈代表該水平下的一個分類，小圓圈直徑大小與相對豐度大小呈正比州刽。著色原則：無顯著差異的物種統(tǒng)一著色為黃色空执，差異物種Biomarker跟隨組進(jìn)行著色，紅色節(jié)點表示在紅色組別中起到重要作用的微生物類群穗椅，綠色節(jié)點表示在綠色組別中起到重要作用的微生物類群辨绊，若圖中某一組缺失，則表明此組中并無差異顯著的物種匹表，故此組缺失门坷。圖中英文字母表示的物種名稱在右側(cè)圖例中進(jìn)行展示。

step1: install lefse through conda

conda create -n lefse -c biobakery lefse -y 
conda activate lefse 
which format_input.py

step2: collapse the table.gza to the L6 level

qiime taxa collapse \
    --i-table result/final_table.qza \
    --o-collapsed-table collapse/collapse.table.qza \
    --p-level 6 \
    --i-taxonomy result/final_taxonomy_sliva.qza

step3: calculate relative-frequency for the collapsed table (relative abundance instead of counts)

qiime feature-table relative-frequency \
    --i-table collapse/collapse.table.qza \
    --o-relative-frequency-table collapse/collapse.frequency.table.qza \
    --output-dir collapse/ 

step4: export biom file

qiime tools export \
    --input-path collapse/collapse.frequency.table.qza \
    --output-path collapse/

step5: convert biom to text file

biom convert \
    -i collapse/feature-table.biom \
    -o collapse/collapse.frequency.table.tsv \
    --header-key "taxonomy" \
    --to-tsv

step6: filter tax

sed 's/;/\|/g' collapse/collapse.frequency.table.tsv | \
    awk '{split($1, a, "|");if( a[6] != "__"){print $0}}' | \
    #sed 's/d\_\_Bacteria|//g' | \
    grep -vE "g__uncultured|d__Archaea|p__WPS-2|p__SAR324_clade|Constructed" | \
    sed 's/#OTU ID/Group/g;s/taxonomy//g' > collapse/collapse.frequency.table.lefse.tsv

step7: run lefse

conda activate lefse
# convert text file into lefse.input file 
format_input.py \
    collapse/collapse.frequency.table.lefse.tsv \
    result/collapse.frequency.table.lefse.in \
    -c 1 \
    -m f \
    -o 100000
# run lefse
run_lefse.py \
    result/collapse.frequency.table.lefse.in \
    result/collapse.frequency.table.lefse.res 
# select significant result Lefse
grep -E "HTN|Normal" \
        result/collapse.frequency.table.lefse.res \
        > result/collapse.frequency.table.lefse_signif.res
# plot lda 
plot_res.py \
    result/collapse.frequency.table.lefse_signif.res \
    result/lefse_final_lda.pdf \
    --format pdf \
    --autoscale 0
# plot cladogram 
plot_cladogram.py \
    result/collapse.frequency.table.lefse_signif.res \
    result/lefse_total_clado.pdf \
    --format pdf

Functional prediction: picrust2

Picrust是Phylogenetic Investigationof Communities by Reconstruction of Unobserved States的簡稱桑孩，是一款基于16s rRNA基因序列預(yù)測微生物群落功能的軟件拜鹤。

其原理：

（1）基因內(nèi)容預(yù)測（gene content inference）。該步先對Greengenes數(shù)據(jù)庫的“closed reference”序列劃分OTU后構(gòu)建進(jìn)化樹流椒，通過祖先狀態(tài)重構(gòu)（Ancestralstate reconstruction）算法并結(jié)合IMG/M數(shù)據(jù)庫敏簿，預(yù)測出樹中未進(jìn)行全基因組測序OTU的基因組信息。

（2）宏基因組預(yù)測（metagenome inference）。將16SrDNA測序結(jié)果與Greengenes數(shù)據(jù)庫進(jìn)行比對惯裕，挑選出與“closed reference”數(shù)據(jù)庫相似性高的（默認(rèn)為≥97%）OTU温数；根據(jù)OTU對應(yīng)基因組中16SrDNA的拷貝數(shù)信息，將每個OTU對應(yīng)序列數(shù)除以其16S拷貝數(shù)來進(jìn)行標(biāo)準(zhǔn)化蜻势；最后撑刺，將標(biāo)準(zhǔn)化的數(shù)據(jù)乘以其對應(yīng)的基因組中基因含量從而實現(xiàn)宏基因組預(yù)測的目的。獲得的預(yù)測結(jié)果可以通過KEGG Orthology握玛、COGs或Pfams等對基因家族進(jìn)行分類够傍。

qiime2-2020.8版本暫時無法安裝q2-picrust插件，因此使用picurst2軟件做微生物功能預(yù)測分析挠铲。

# install picrust2
conda create -n picrust2 -c bioconda -c conda-forge picrust2=2.3.0_b -y 
# export representative sequences
conda activate qiime2-2020.8
qiime tools export --input-path result/final_rep_seqs.qza --output-path ./ 
conda deactivate 
# run picrust2
conda activate picrust2
picrust2_pipeline.py -s dna-sequences.fasta -i feature-table.biom -o picrust2_out_pipeline -p 30
conda deactivate

The key output files are:

• EC_metagenome_out - Folder containing unstratified EC number metagenome predictions (pred_metagenome_unstrat.tsv.gz), sequence table normalized by predicted 16S copy number abundances (seqtab_norm.tsv.gz), and the per-sample NSTI values weighted by the abundance of each ASV (weighted_nsti.tsv.gz).
• KO_metagenome_out - As EC_metagenome_out above, but for KO metagenomes.
• pathways_out - Folder containing predicted pathway abundances and coverages per-sample, based on predicted EC number abundances.

Reference

1. qiime2 docs
2. fastqc tutorial
3. usearch tutorial
4. vsearch tutorial
5. import data in qiime2
6. 擴增子分析軟件qiime2必知必會
7. 擴增子數(shù)據(jù)庫整理
8. Greengene數(shù)據(jù)庫整理
9. SILVA 數(shù)據(jù)庫整理
10. 差異檢驗
11. lefse after qiime2
12. lefse scripts github
13. picrust2安裝及對16s數(shù)據(jù)進(jìn)行功能預(yù)測
14. picrust2 wiki