# Qiime2環境安裝

**1. 下載QIIME 2 conda配置檔.yml:**

wget https://data.qiime2.org/distro/core/qiime2-2023.5-py38-linux-conda.yml

#

**2. 編輯配置檔加入nextflow與krona**

- nextflow=23.04.1
- krona=2.8.1

#

**3.啟動conda並建立環境**

source /storage_1/KF07453/miniconda3/bin/activate
conda env create -n qiime2 --file qiime2-2023.5-py38-linux-conda.yml

#

**4. 顯示所有環境並啟動環境**

conda env list
conda activate qiime2

#

**5.環境設定**
將新建立的qiime2環境的binary執行檔路徑寫入環境變數中
先進入python環境--->python

import os
Add_Binarry_Path='~/.conda/envs/HMP_project/bin/'
os.environ['PATH']=Add_Binarry_Path+':'+os.environ['PATH']
exit()

#
**6.建立專案目錄**
project_folder = "mouse_tutorial"

mkdir -p ./mouse_tutorial/seqs
mkdir -p ./mouse_tutorial/qza
mkdir -p ./mouse_tutorial/qzv
cd mouse_tutorial

# QIIME2分析-資料
Qiime2為用來分析微生物次世代基因定序資料的開源軟體

**1. 測試資料說明**
本次測試資料源自於Qiime2官網文件的Parkinson's Mouse Tutorial教學

使用之資料源自於Sampson et al., 2016的研究
(https://www.cell.com/cell/fulltext/S0092-8674(16)31590-2?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0092867416315902%3Fshowall%3Dtrue)
此實驗證實了糞便之微生物菌相對Parkinson’s Disease (PD)症狀發展有所影響

本次測試資料來自其中部分實驗資料，探討移植之菌相來源對於小鼠之微生物菌相分布之影響，而菌相分別源自健康對照組及PD患者之糞便微生物

**2.下載測試資料**

wget \
-O "demultiplexed_seqs.zip" \
"https://data.qiime2.org/2023.5/tutorials/pd-mice/demultiplexed_seqs.zip"

**3.將測試資料解壓縮並移動到seq資料夾中**

unzip demultiplexed_seqs.zip
mv demultiplexed_seqs/* seqs

**4.建立Manifest樣本清單檔案(manifest.tsv)**
Manifest.tsv目的是告訴Qiime2 .fastq.gz的路徑在哪裡，以及該檔案所對應的樣本名稱

可使用manifest.py生成manifest.tsv，或使用excel製作
結果會顯示:
manifest.tsv already generate in /storage_1/KF07453/mouse_tutorial/

查看manifest.tsv的前10行

head -n 10 manifest.tsv


**5.下載metadata詮釋資料**
metadata記錄實驗各樣本之不同特徵
若需要更改metadata.tsv，可用excel進行修改

wget \
-O "metadata.tsv" \
"https://data.qiime2.org/2023.5/tutorials/pd-mice/sample_metadata.tsv" 

查看metadata.tsv

head -n 10 metadata.tsv


columns儲存樣本名稱(sample_name)，barcode，實驗小鼠編號(mouse_id)，基因型(genotype)，籠子編號(cage_id)，菌相來源(donor)，人類健康狀態(donor_status)，微生物菌相移植天數(days_post_transplant)，小鼠基因型與人類健康狀態(genotype_and_donor_status)。
Row 2紀錄每個特徵的屬性，樣本名稱(column 1)會用#q2:types表示，其餘的依照類別型(categorical)或數值型(numeric)進行分類。接著每個row代表一個樣本
#

# QIIME2 分析
**1.輸入原始資料，儲存成qza格式**

qza: Qiime2中儲存資料的格式，任何分析時input皆為qza檔案

qzv: 與qza對應的格式，供視覺化的檔案，可上傳至qiime2 view上查看結果

qiime tools import \
--type "SampleData[SequencesWithQuality]" \
--input-format SingleEndFastqManifestPhred33V2 \
--input-path ./manifest.tsv \
--output-path ./qza/demux_seqs.qza

'''
--type "SampleData[PairedEndSequencesWithQuality]" \
--input-format PairedEndFastqManifestPhred33V2 \
--input-path ./manifest.tsv \
--output-path ./qza/demux_seqs.qza
'''


* qiime tools import: 用此命令將資料輸入
* --type: 資料的形式，
單端序列格式:SampleData[SequencesWithQuality]
成對定序格式:SampleData[PairedEndSequencesWithQuality]
* --input-format:Fastq格式中序列品質的編碼方式，主流為Phred33V2，其次還有Phred64V2
若資料為成對定序並以Phred33V2進行編碼，則指定--input-format PairedEndFastqManifestPhred33V2
* --input-path:manifest之檔案路徑
* --output-path:qza檔案輸出的位置，此範例設定輸出在當前目錄下，並命名為demux_seqs.qza
#

**2.檢視樣本的定序深度(sequencing depth)與序列的品質(quality score)**

qiime demux summarize \
--i-data ./qza/demux_seqs.qza \
--p-n 20000 \
--o-visualization ./qzv/demux_seqs.qzv


* --i-data: demux_seqs.qza的檔案位置
* --p-n: 從所有序列中隨機選擇的序列數量，預設為10000，選中的序列用以生成序列品質分布圖(quality plot)
* --o-visualization: qzv檔案的輸出位置，命名為demux_seqs.qzv

將demux_seqs.qzv檔案上傳到qiime2 view(https://view.qiime2.org/)進行視覺化

從quality plot可判斷序列的整體品質分布，並決定是否要進行修剪(trimmed)，以及從哪個位置進行修剪(position)
將游標移動到box，會顯示該位置的seven-number summary
通常定序剛開始與結束的兩端品質都會比較差，若序列頭尾品質低於20可考慮切除



Q = -10log10p（p表示測序的錯誤率，Q表示鹼基質量分數）

Q Value,正確率,錯誤率
10, 90%, 10%
20, 99%, 1%
30, 99.9%, 0.1%

#

**3.序列品質控制(DADA2)**
DADA2套件: 可對序列進行修剪(trimming)、過濾(filters)、降噪(denoises)、合併(merge)、去除重疊(dereplicates)

因上圖quality plot顯示在任何位置，Q value幾乎大於20，定序品質良好，因此保留完整序列(150 bases)

qiime dada2 denoise-single \
--i-demultiplexed-seqs ./qza/demux_seqs.qza \
--p-trim-left 0 \
--p-trunc-len 150 \
--p-n-threads 0 \
--o-denoising-stats ./qza/dada2_stats.qza \
--o-table ./qza/dada2_table.qza \
--o-representative-sequences ./qza/dada2_rep_set.qza


*qiime dada2 denoise-single: 定序品質管控指令，輸出3個qza檔案
*--i-demultiplexed-seqs: 輸入之序列路徑
*--p-trim-left: 序列起始端切多少，預設是0
*--p-trunc-len: 序列剪切後的總長度，若單一序列低於此設定長度則捨棄這一序列
*--p-n-threads: 使用的CPU數量，設定0則使用所有的CPU(預設為1)
*--o-denoising-stats: 輸出檔案的位置，使用dada2各個階段控制序列的簡單統計(statistics)。檢查各個控管階段有無不尋常的篩選掉大量序列，以評估是否剪切過多
*--o-table: dada2處理後之特徵表(feature table)，包含了每個樣本中的微生物特徵（如ASVs或OTUs）的計數信息，以Feature ID表示描述這些特徵在不同樣本中的存在。特徵表對於後續的微生物群落分析非常重要，它可以被用來比較不同樣本之間的微生物組成，進行多樣性分析
*--o-representative-sequences: 各個Feature ID所對應的代表序列(representative sequences)

*qiime dada2 denoise-paired: 雙尾定序品質控管指令
*--p-trunc-len-f ...: forward-序列剪切後的總長度，若單一序列低於此設定長度則捨棄這一序列
*--p-trunc-len-r ...: reverse-序列剪切後的總長度，若單一序列低於此設定長度則捨棄這一序列

#================================================================

經序列品質管控後，會將相同的序列分成1組ASV（amplicon sequence variant)，並從各ASVs中選出一條代表該群的序列(Representative sequences)，並給予一串標號(Feature IDs)

ASV(擴增子序列變異) vs OTU(操作分類單元)

OTU: 跟據序列的相似性，以一定的相似度把序列歸集到一起。一般用的相似度是 97% 或 99%，產生出來的OTU稱為OTU97或 OTU99。
缺點: 敏感度不足，一個 OTU 中可能包含有不同的實際物種；再者，如果序列因為定序過程產生錯誤，會造成不可信賴的聚類結果

ASV: 針對定序時的錯誤來做修正。以目前最主流的 DADA2 演算法來說，就是根據 Illumina 定序的錯誤模型，藉由演算法與非監督式學習方法 selfConsist 來去除因為定序所造成的序列錯誤，而產出 ASV。因此相較於 OTU，ASV 的區分精細度就比較高。

#================================================================

#

**4.qza轉成qzv視覺化**
(1)表格化降噪統計結果

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


表格化統計結果，input(輸入數據)；filtered(過濾後的數據)；denoised(降噪後的數據)；non-chimeric(去除嵌合體後的數據)

(2)特徵表結果
查看每個樣本品質控管後的特徵量，以及每個特徵頻率和其在樣本中出現的次數

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

!! 若有出現以下Error，代表**樣本名稱與metadata.tsv不同**，主要為上傳到qiime2時，sample-id多了_和後面的數字

解決方法：可將qzv資料夾中的dada2_stats.qzv上傳到qiime2 view，將視覺化的檔案下載(download metadata TSV file)，並用excel打開，複製sample-id並取代原本的metadata.tsv的sample-name


最後再執行一次指令即可


(3)特徵ID的代表序列

qiime feature-table tabulate-seqs \
--i-data qza/dada2_rep_set.qza \
--o-visualization qzv/dada2_rep_set.qzv


**5.物種分類(Taxonomy assignment)**
使用機器學習分類器進行物種分類

* 分類器(classifier): 用以分類菌種之機器學習模型
* 微生物序列資料庫(microbiota database): 用以訓練分類器之資料來源，常見的微生物資料庫，分別為Greengenes, SILVA和RDP，下圖擷取自https://ithelp.ithome.com.tw/articles/10298358
* 已訓練的分類器:有訓練好的分類器模型可直接使用(https://docs.qiime2.org/2023.5/data-resources/#taxonomy-classifiers-for-use-with-q2-feature-classifier)


(1)下載物種分類器

* 本次分析的序列為16S V4，將使用qiime2提供之Silva 138 99% OTUs from 515F/806R region of sequences分類器，可透過`wget`下載分類器，並命名為`silva-v4-nb-classifier.qza`

wget \
-O "qza/silva-V4-nb-classifier.qza" \
"https://data.qiime2.org/2023.5/common/silva-138-99-515-806-nb-classifier.qza"


(2)進行物種分類

qiime feature-classifier classify-sklearn \
--i-classifier qza/silva-V4-nb-classifier.qza \
--i-reads qza/dada2_rep_set.qza \
--o-classification qza/taxonomy.qza

* qiime feature-classifier classify-sklearn: 透過分類器對representive sequences進行分類之指令
* --i-classifier: 分類器位置
* --i-reads: representive sequences之檔案位置
* --o-classification: 輸出結果之taxonomy.qza

(3)將分類結果視覺化

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



* Taxon 呈現各階層的分類結果 (界d、門p 綱c、目o、科f、屬g、種s)
* 若是出現許多顯示 Unassigned 或多數無法到 G (Genus) / S (Species)， 就要考慮分類器是否合適。由於使用16S V4定序，片段較短，許多分類結果只能達到G階層。

#
**6.分類柱狀圖**

qiime taxa barplot \
--i-table qza/dada2_table.qza \
--i-taxonomy qza/taxonomy.qza \
--m-metadata-file metadata.tsv \
--o-visualization qzv/taxa-bar-plots.qzv


* Download: 可選擇資料之下載方式，圖片(柱狀圖bar or 資料標籤legand)或是csv(儲存不同樣本之metadata及在特定分類層級下，各微生物出現頻率)
* Bar Width: 可調整每個bar的寬度
* Taxomic Level: 可選擇繪製之分類層級(界1 門2 綱3 目4 科5 屬6 種7)
* Color Palette: 可選擇bar繪製顏色
* Sort Sample By: 可選擇樣本依照metadata分類之排序依據


不同樣本之分類柱狀圖，顯示不同樣本之微生物相對頻率(Relative Frequency)，也可以想成是微生物於樣本中之相對豐度(Relative Abundence)。顏色代表微生物分類，分類層級將依照使用者自行選擇進行繪製(這邊是選擇level 4 -->目)，游標移動至柱狀圖上將顯示詳細微生物在樣本中佔比於圖片上方

#
**7.互動式動態分類圓餅圖(krona)**

使用krona套件，需另外透過conda安裝korna

`conda install -c bioconda krona`

並透過pip安裝q2-krona
`pip install git+https://github.com/kaanb93/q2-krona.git`

qiime krona collapse-and-plot \
--i-table qza/dada2_table.qza \
--i-taxonomy qza/taxonomy.qza \
--o-krona-plot qzv/krona.qzv

* qiime krona collapse-and-plot: 繪製動態分類圓餅圖指令
* --i-table:輸入之feature table qza檔案位置
* --i-traxonomy: 輸入之taxonomy.qza檔案位置
* --o-krona-plot: 輸出之krona.qzv位置


圓餅圖表示單一樣本中微生物在不同階層之分布，可以透過游標按住並滑動來觀察特定分類階層下的詳細分類
最內層代表界層級之分類，最外圍則為種層級之分類(根據taxonamy之分類解析度，最外層不一定是種)


**8.繪製熱圖**

qiime taxa collapse \
--i-table qza/dada2_table.qza \
--i-taxonomy qza/taxonomy.qza \
--p-level 4 \
--o-collapsed-table qza/col_table_4.qza

* qiime feature-table heatmap: 繪製熱圖之指令
* --o-visualization qzv: 輸出之heatmap位置
* --p-color-scheme: 指定heatmap的配色(https://matplotlib.org/stable/users/explain/colors/colormaps.html)
* --p-cluster: 分群，可以選擇`features`、`samples`、`both`、`none`，根據序列/樣本/兩者/都不要


進行聚類(clustering)之結果


這張heatmap是依照samples進行cluster，並且依照level4進行collaspe的結果
是根據樣本進行cluster的，因此若是樣本有依不同處理進行分組的話，就能快速地看出處理組之間的菌像分布較相似
若是依照features進行cluster圖左的黑線則會依照feature相似性顯示在heatmap上方

#

**9.繪製親緣樹(phylogenetic tree)**

分為rooted tree以及unrooted tree

rooted tree: 物種起源皆自同一個點，可看出演化的方向性，適合應用在觀察演化過程中不同物種之演化關係

unrooted tree: 定義演化的方向與路徑，著重在分類群之間的關係

繪製親緣數除了供視覺化參考外，後續多樣性分析的部分指標會需要用rooted tree來考慮菌種間的親緣關係(如Faith's Phylogenetic Diversity和UniFrac distance)

樹的製作主要可分為兩種方法:
* **fragment-insertion**: 將representative sequences依照序列插入已知的參考親緣樹中，需較多的計算和時間成本
* **align-to-tree**: 將序列進行alignment，接著透過遮罩(masks)，去除高度變異的位置，最後繪製親緣樹

qiime2官方教學文件解釋，fragment-insertion在對短片段Illumina序列進行親緣樹時可能會比傳統的基於序列比對的方法表現更好，儘管會需要更多的計算資源，因此本範例使用fragment-insertion方法繪製親緣樹

(1)下載參考樹

wget \
-O "qza/sepp-refs-gg-13-8.qza" \
"https://data.qiime2.org/2023.5/common/sepp-refs-gg-13-8.qza"


(2)進行親緣樹比對

qiime fragment-insertion sepp \
--i-representative-sequences ./qza/dada2_rep_set.qza \
--i-reference-database qza/sepp-refs-gg-13-8.qza \
--o-tree ./qza/rooted_tree.qza \
--o-placements ./qza/tree_placements.qza \
--p-threads 80

* qiime fragment-insertion sepp: 透過sepp方法進行fragment-insertion之指令
* --i-representative-sequences: 輸入之representative sequences qza檔案位置
* --i-reference-database: 輸入之參考樹位置
* --o-tree: 輸出之tree位置及名稱
* --o-placements: 輸出之插入結果之文件
* --p-thread: 使用之cpu核心數量，可透過`nproc`查詢最大CPU可使用量


**10.多樣性分析**
多樣性分析分為alpha及beta多樣性

* alpha多樣性: 反應**單一樣本(或1組樣本)內之多樣性**，不涉及個體之間的比較
--> 物種豐富度(richness)、物種均勻度(evenness)

* beta多樣性: 組間比較(兩組/多組間菌種分布差異程度)
以下3種多樣性的指標的視覺化結果可分為distance和pcoa(以emperor)

--> bray curtis: 考量物種組成也採計各物種的數量(abundance)-->物種共有的個體數與兩群集規模平均的比值
--> jaccard: 只比較兩群及共有物種的數量(richness)-->交集/聯集，數值越高，樣本相似越高
--> unifrac: 特定群集在系統發生樹上獨佔節點的比例，是評估群集中物種親緣關係的多樣性指標，基於物種的進化關係，判斷樣本間的相似程度
又分為weighted和unweighted，對應是否採計族群規模進行加權

#

**11.Alpha多樣性之稀疏曲線(rarefaction curve)繪製**

稀疏曲線的繪製: 需有各樣本ASVs資訊(feature table) + 各ASVs間的親緣關係(phylogenetic tree)

(1)對各樣本進行標準化
經標準化後，才能比較不同樣本之間多樣之差異

qiime diversity alpha-rarefaction \
--i-table ./qza/dada2_table.qza \
--i-phylogeny qza/rooted_tree.qza \
--m-metadata-file ./metadata.tsv \
--o-visualization ./qzv/alpha_rarefaction_curves.qzv \
--p-min-depth 1 \
--p-max-depth 4996


* qiime diversity alpha-rarafaction: 繪製alpha rarefaction之指令
* --i-table: 輸入之feature table qza檔案之位置
* --i-phylogeny: 輸入之rooted tree qza檔案位置，要輸入親緣樹才能計算Faith's phylogenetic diversity
* --i-metadata-file: 輸入之metadata之位置
* --o-visualization: 輸出之qzv結果位置
* --p-min-depth: 繪製稀疏曲線之最小測序深度，預設為1
* --p-max-depth: 繪製稀疏曲線之最大測序深度，可查看feature table qzv設定樣本之最大測序深度(此範例為4996)
#

* Metrix: 可選擇alpha多樣性指標
樣本中觀察到物種數量observed_features在特定採樣深度下，觀察到的ASVs數量
* Shannon多樣性指標和系統發育多樣性faith_pd: 皆為多樣性指數，計算樣本之物種豐富度(richness)和均勻度(evenness)，差別在於faith_pd計算時有考慮到物種間親緣關係(phylogenetic tree)
* Sample Metadata column: 繪製稀疏曲線樣本分組方式


* X軸: 採樣深度(sampling depth)
* y軸: 多樣性指標
* 可透過稀疏曲線確認**樣本中的序列數是否足以代表整個樣本**，取樣深度代表每一輪無放回抽樣時所取出的菌種觀察值數量，隨採樣深度提高，Y軸之多樣性指標也趨於飽和(斜率接近於0)，即可推論樣本中的序列數足以代表整個樣本


* 上圖顯示每個取樣深度下每個分組的樣本數，有助於確定在特定取樣深度下的樣本損失
* 深度的選擇在於最小化序列損失的同時，最大化保留用於多樣性分析的序列，本範例將使勇每個樣本2000個序列的取樣深度進行分析


* qiime diversity core-metrics-phylogenetic: 繪製alpha多樣性之指令
* --i-table: 輸入之feature table qza檔案之位置
* --i-phylogeny: 輸入之rooted tree qza檔案位置，有輸入親緣樹才能計算Faith's phylogenetic diversity
* --m-metadata-file: 輸入之metadata之位置
* --p-sampling-depth: 取樣深度
* --output-dir: 將所有結果輸出成一個資料夾

**12.Alpha多樣性**
alpha多樣性使用faith_pd，jaccard和bray curtitis的視覺化qzv檔案已存在core-metrics-results資料夾中

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


* qiime diversity alpha-group-significance: 視覺化alpha多樣性之指令
* --i-alpha-diversity: 輸入之alpha多樣性檔案位置
* --m-metadata-file: 輸入之metadata之位置
* --o-visualization: 輸出之檔案位置


* 圖片上方Column為分組依據
* 顯著差異: 探討不同條件下分組之組間差異，並且進行成對比較(pairwise)檢測差異是否顯著

將core-metrics-results資料中之bray_curtis_emperor.qzv透過emperor套件視覺化

#

**13.Beta多樣性視覺化**

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


* qiime diversity beta-group-significance: 視覺化beta多樣性之指令
* --i-distance-matrix: 輸入之beta多樣性(for distance)檔案位置
* --m-metadata-file: 輸入之metadata位置
* --m-metadata-column: 樣本分組依據
* --p-pairwise: 進行成對比較(pairwise)，預設只有ANOVA
* --o-visualization: 輸出檔案之位置

分組(donor type)之間(PERMANOVA, Permutational analysis of variance)分析結果
上圖按照genotype_and_donor_status分組的4組小鼠之間微生物beta多樣性(distance)存在統計上之顯著差異

* PERANOVA: 是一種比較兩分組之間的距離，為一種無母數統計的分析方法
* Number of permutations: 無母數統計中，對數據進行排序的方式的數量，用於計算p值
#


* 樣本共分為4組，因此會產生4張圖，這邊只顯示其中2個結果，2張圖表示4個組別分別與susceptible and Healthy以及susceptible and PD兩組間的距離盒鬚圖，數值越大代表兩組之間差異越大
* 確定PERMANOVA分析組間的距離存在顯著差異後，進行成對比較可看出哪兩組間有顯著差異，可看出四組之間距離皆存在顯著差異
* 可將--m-metadata-column改成--m-metadata-column donor或--m-metadata-column genotype查看更詳細的分析結果

**Emperor**
Emperor是強大的互動式PCoA(主座標分析)視覺化工具，功能有很多，可選擇樣本的分組方式(透過**顏色**或是**形狀**)，也可選擇顯示特定分類的數據(Visibility)，並可選擇顯示的主座標(Axes)


* 主成分分析(PCA)與主座標分析(PCoA)皆為一種資料降維方法，對多為數據進行降維，同時保持數據中雙方差異貢獻最大的特徵，若樣本的群落組成越相似，則PCA/PCoA圖中的距離則越接近，因此群落結構相似度高的樣本會聚集在一起
* PCA: 在高維數據中找到最能解釋變異性的主成分(設法保留數據裡的變異(variance)，讓點的位置盡量不要變動)
* PCoA: 主要使用距離矩陣，在多維空間中呈現樣本之間的距離關係(讓原本數據各點間距哩，跟投影裡各點間距離之間的相關最高)


#

**14.差異豐富度分析(LDA Effect size, LEFse)**

LEFse可用來檢測不同樣本組特徵(如ASV、菌種)之間的相對豐富度是否存在顯著差異

(1)過濾掉低豐富度(abundance)與低出現頻率(prevalence)的ASVs，可提供更好的分辨率

qiime feature-table filter-features \
--i-table ./qza/dada2_table.qza \
--p-min-frequency 50 \
--p-min-samples 4 \
--o-filtered-table ./qza/dada2_table_filtered.qza


* qiime feature-table filter-features: 過濾特定AVs的指令
* --i-table: 輸入feature table之位置
* --p-mjin-frequency: ASV在所有樣本中出現的總頻率低於這個值，則去除掉此ASV
* --p-min-samples: ASV在樣本中出現的次數，低於這個值，則去除掉此ASV
* --o-filtered-table: 輸出篩選後的filtered table
#
* --m-metadata-file: 輸入metadata之路徑，本範例未使用
* --p-where: 若有輸入metadata，則可設定此參數透過metadata分組，篩選特定分組的結果，如donor=hc

(2)安裝dokdo套件

先透過git取得GitHub的dokdo檔案，再透過pip進行安裝，而不能像krona或nextflow直接先寫在yml配置檔裡面，因為pip還存在著另一個dokdo套件，兩者是不一樣的東西，以避免與我們要的dokdo搞混

git clone https://github.com/sbslee/dokdo
cd dokdo
pip install .
cd ..


(3)產生LEFse

dokdo prepare-lefse \
--table-file ./qza/dada2_table_filtered.qza \
--taxonomy-file ./qza/taxonomy.qza \
--metadata-file ./metadata.tsv \
--class-col donor \
--output-file ./lefse_input_table.tsv


* dokdo prepare-lefse: 輸出LEFse之指令
* --table-file: 輸入feature table之位置，本範例使用的是filtered table
* --taxonomy-file: 輸入taxonomy.qza之位置
* --m-metadata-file: metadata之檔案路徑

(4)透過網頁版[Galaxy with LEFse](https://galaxy.biobakery.org/)進行視覺化

(i)上傳檔案，點選左上角的上傳符號

(ii)按下方的choose local file-->選擇lefse_input_table.tsv上傳


(iii)按左方的工具列的LEfSe，從A(按下A-->上傳檔案-->執行)到C


#

**15.微生物基因功能預測PICRUSt2**

PICRUSt2(Phylogenetic investigation of communities by reconstruction of unobserved states)為一套根據標記基因的序列(如16S rRNA)去預測微生物基因功能的工具

PICRUSt2會將各個ASV序列與其參考資料中的序列進行比對，推斷每個ASV所含的基因家族(Gene family)，從而預測功能


(1)將dada2_table.qza透過qiime2輸出成feature-table.biom

Biom(Biological observation matrix)檔案是專為生物資訊分析的檔案，內含ASVs在各樣本中的數量(abundance)資訊

qiime tools export \
--input-path ./qza/dada2_table_filtered.qza \
--output-path ./


* qiime tools export: 用於將qza匯出成不同格式
* --input-path: 輸入feature table的位置，本範例使用的是filtered table
* --output-path: 輸出之BIOM位置

(2)檢視輸出的BIOM檔案

`biom head -i feature-table.biom`


(3)將dada2_rep_set.qza轉成fasta，內含feature ID與原始序列

qiime tools export \
--input-path qza/dada2_rep_set.qza \
--output-path ./


(4)檢視輸出的fasta檔案

`head dna-sequences.fasta`


#
**16.安裝與啟動mamba及PICRUSt2 (v2.5.0)**

官方建議設備: 含有16GB以上的RAM
安裝: 使用mamba而非conda，安裝速度會比conda還快，且使用conda安裝picrust2會產生套件之間的版本衝突

**(1)安裝mamba與picrust2環境**

conda deactivate
conda install mamba -n base -c conda-forge

mamba create -n picrust2 -c bioconda -c conda-forge picrust2=2.5.0

conda activate picrust2




**(2)將安裝完成之picrust2寫入環境變數**

python
import os
Add_Binarry_Path='~/.conda/envs/picrust2/bin/'
os.environ['PATH']=Add_Binarry_Path+':'+os.environ['PATH']
exit()

**(3)執行Picrust2分析**
執行picrust2_pipeline.py

picrust2_pipeline.py \
-s dna-sequences.fasta \
-i feature-table.biom \
-o picrust2_out_pipeline \
-p 80


* -s: 輸入之sequence.fasta位置
* -i: 輸入之feature-table.biom位置
* -o: 輸出檔案之資料夾名稱
* -p: CPUs數量

若想要重跑Picrust2，輸出結果至同樣命名的folder，則要刪掉原本的folder才能重跑

**(4)進入picrust2_out_pipeline資料夾**

`cd picrust2_out_pipeline`

進入picrust2_out_pipeline資料夾，可以看到4個資料夾和4個檔案


* EC_metagenome_out: 使用MetaCyc代謝路徑資料庫的酵素編號(enzyme commission, EC numbers)豐富度預測結果、每個樣本的NSTI權重(weighted_nsti.tsv)、預測16s校正的特徵表
* intermediate: Folder containing intermediate MinPath files and files used for sequence placement pipeline
* KO_metagenome_out: 使用KEGG生物路徑資料庫的KO基因豐富度預測結果
* pathways_out: 根據EC numbers豐富度預測結果所推論出的MetaCyc代謝路徑豐富度
* EC_predicted.tsv.gz: 預測每個ASV中EC的數量
* KO_predicted.tsv.gz:預測每個ASV中對應到的KEGG pathway數量
* marker_predicted_and nsti.tsv.gz: 預測每個ASV中16s與NSTI的數量
* out.tre: 參考序列的樹

#

-->Picrust2預設使用MetaCyc資料庫進行pathway預測，因MetaCyc亦是微生物總體基因體定序(MGS, Shotgun metagenomics)基因功能分析常見的資料庫

-->Picrust2可藉由菌相資料獲得Enzyme Commission(EC number)、KEGG、MetaCyc路徑豐富度資料，以此預測菌群中的酵素及代謝路徑豐富度

-->Picrust2 Pathway的豐富度計算以預測EC number 的豐富度為主，並基於EC number比對MetaCyc資料庫進而計算出 Pathway的豐富度以降低產出假陽性 Pathways 的可能性

-->得到功能預測結果後，常會進行組間的豐富度比較，官方建議可採用具使用者介面的STAMP或是ALEDx2方法進行統計檢驗

**(5)輸出加上註釋的tsv檔案**

-->分析後主要透過三個壓縮檔解壓縮後的tsv視覺化，分別為
`EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz`
`KO_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz`
`pathways_out/path_abun_unstrat_descrip.tsv.gz`

然而Picrust2輸出時僅含有酵素及代謝路徑編號，為方便人類判讀，還會使用[Picrust2 Add description](https://github.com/picrust/picrust2/wiki/Add-descriptions)功能，替編號加入註釋(EC number pathway)，再加上註釋的tsv忠看到第二個column為description，供人類方便檢視

add_descriptions.py \
-i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz \
-m EC \
-o EC_metagenome_out/EC_pred_metagenome_unstrat_descrip.tsv.gz

!gunzip EC_metagenome_out/EC_pred_metagenome_unstrat_descrip.tsv.gz

add_descriptions.py \
-i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz \
-m KO \
-o KO_metagenome_out/KEGG_pred_metagenome_unstrat_descrip.tsv.gz

!gunzip KO_metagenome_out/KEGG_pred_metagenome_unstrat_descrip.tsv.gz

add_descriptions.py \
-i pathways_out/path_abun_unstrat.tsv.gz \
-m METACYC \
-o pathways_out/path_abun_unstrat_descrip.tsv.gz

!gunzip pathways_out/path_abun_unstrat_descrip.tsv.gz


* add_descriptions.py: 於tsv加上註釋
* -i: 輸入之tsv.gz檔案位置
* -m: 使用的註釋
* -o: 輸出檔案之名稱
* gunzip: 解壓縮.gz

#

**17.STAMP分析**
[STAMP](https://picrust.github.io/picrust/tutorials/stamp_tutorial.html)可對tsv檔案進行視覺化，以查看預測菌群中的酵素及代謝路徑豐富度之結果

(1)於本機端下載[STAMP](https://beikolab.cs.dal.ca/software/STAMP)

(2)將picrust2輸出的tsv檔案下載到本機端，並以Excel開啟

(3)刪除EC_pred_metagenome_unstrat_descrip.tsv與KEGG_pred_metagenome_unstrat_descrip.tsv中的function、path_abun_unstrat_descrip.tsv中的pathway欄(否則輸入STAMP時會發生Error)

(4)將metadata.tsv下載到本機端

(5)打開Stamp，輸入`ctrl+O`載入檔案
-->Profile file: KO/EC/Pathway之tsv
-->Group metadata file: metadata.tsv


-->1. 選擇分組依據
-->2. 選擇Multiple groups, two groups, two samples進行比較，並選擇統計方式
-->3. 選擇製作之圖表類型
-->4. 選擇是否過濾部分樣本

#
範例:

Multiple groups / bar plot / 1,4-dihydroxy-2-naphthoate polyprenyltransferase。顯示不同樣本中1,4-dihydroxy-2-naphthoate之豐富度，並以顏色區分不同組

ANOVA統計之P值顯示在圖片右上角。


#

Multiple groups / box plot / 1,4-dihydroxy-2-naphthoate。同上，但是是將組別內所有樣本集合起來以box plot顯示。


#

Multiple groups / post hoc / 1,4-dihydroxy-2-naphthoate
當ANOVA分析結果顯示多組間存在顯著差異，再透過post hoc查看哪些組別存在差異


#

Two groups/ Extended error bar，顯示2組別間具有顯著差異的特徵