內(nèi)部培訓(xùn)資料（2022.10.29）
QIIME 2分析實(shí)例--人類腸癌不同部位微生物多樣性分析
“Moving Pictures” tutorial
熊金波實(shí)驗(yàn)室出品 ??????內(nèi)部交流培訓(xùn)使用
寧波大學(xué)海洋學(xué)院 Larry 陸
本節(jié)1.6萬字堤器，14張圖。閱讀時(shí)間大約50分鐘。
在本次分享中干旧，我們將使用QIIME2軟件進(jìn)行微生物數(shù)據(jù)分析只厘。在五個(gè)時(shí)間點(diǎn)對(duì)來自兩個(gè)人四個(gè)身體部位的微生物組樣本進(jìn)行分析溉躲，第一個(gè)時(shí)間點(diǎn)取樣之后使用了抗生素處理政己∈荒溃基于這批數(shù)據(jù)的研究可以參考2011年發(fā)表在Genome Biology的《Moving pictures of the human microbiome》（10分+的文章哦）杯缺。本次分享將會(huì)完整完成文章中的所有擴(kuò)增子數(shù)據(jù)分析蒸播。本教程中使用的數(shù)據(jù)基于Illumina HiSeq產(chǎn)出，使用地球微生物組計(jì)劃擴(kuò)增16S rRNA基因高變區(qū)4（V4）測(cè)序的方法萍肆。
安裝
Qiime一直以來被詬病最多的就是非常難以安裝袍榆，目前據(jù)我了解有三種主流安裝方法，根據(jù)環(huán)境任選其一即可，實(shí)際的安裝涉及到服務(wù)器運(yùn)維工作，與我們實(shí)驗(yàn)室日常工作相關(guān)性不大徐绑，所以此處不做過多介紹窟绷。我們205實(shí)驗(yàn)室的的T430服務(wù)器和Biostar服務(wù)器都會(huì)部署好QIIME1代和2代，可以直接使用驯用。如果對(duì)于個(gè)人PC有安裝需求的同學(xué)，可以以后向我了解相關(guān)材料。
啟動(dòng)QIIME2運(yùn)行環(huán)境
我們?cè)诜?wù)器上進(jìn)行生物信息數(shù)據(jù)分析的時(shí)候非常重要的就是有良好的運(yùn)行環(huán)境習(xí)慣赞草，有了良好的環(huán)境管理習(xí)慣可以避免大量報(bào)錯(cuò)和數(shù)據(jù)分析問題。
我們?cè)诿看畏治鲩_始前吆鹤，必須先進(jìn)入工作目錄厨疙，除非你是一個(gè)把什么東西都放在桌面上還很工作更有效率的人??。
# 定義工作目錄變量檀头，方便以后多次使用
wd=~/github/QIIME2ChineseManual/2019.7
mkdir -p $wd
# 進(jìn)入工作目錄轰异，是不是很簡(jiǎn)介，這樣無論你在什么位置就可以快速回到項(xiàng)目文件夾
cd $wd
# 方法1\. 進(jìn)入QIIME 2 conda工作環(huán)境
conda activate qiime2-2019.7
# 這時(shí)我們的命令行前面出現(xiàn) (qiime2-2019.7) 表示成功進(jìn)入工作環(huán)境
# 方法2\. conda版本較老用戶暑始，使用source進(jìn)入QIIME 2
source activate qiime2-2019.7
# 方法3\. 如果是docker安裝的請(qǐng)運(yùn)行如下命令搭独，默認(rèn)加載當(dāng)前目錄至/data目錄
docker run --rm -v $(pwd):/data --name=qiime -it qiime2/core:2019.7
# 創(chuàng)建本節(jié)學(xué)習(xí)目錄
mkdir qiime2-moving-pictures-tutorial
cd qiime2-moving-pictures-tutorial
激活分析環(huán)境
qiime2
創(chuàng)建個(gè)人文件夾
mkdir -p larry
首先獲取分析的數(shù)據(jù)
cp -r raw-data larry
cp TruSeq3-PE.fa larry
首先獲取分析樣本元數(shù)據(jù)信息表格
cp mapping_file.txt larry
正式開始分析數(shù)據(jù)要先對(duì)數(shù)據(jù)進(jìn)行指控流程
激活軟件所在的環(huán)境
cd larry
conda activate /biostack/bioconda
創(chuàng)建質(zhì)控文件夾
mkdir -p data_analysis/trimming/reads
參數(shù)解析 :
使用 [trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) 去除拆分序列中可能可能存在的接頭和切掉低質(zhì)量的堿基， 針對(duì)一些小片段擴(kuò)增產(chǎn)物廊镜，測(cè)序數(shù)據(jù)可能包含接頭序列牙肝。
命令模式: SE指定單端數(shù)據(jù) ，PE指定雙端數(shù)據(jù)，對(duì)于PE模式: 2個(gè)輸入文件（正向和反向reads）和4個(gè)輸出文件（正向配對(duì)配椭、正向未配對(duì)虫溜、反向配對(duì)和反向未配對(duì)）
threads： 設(shè)定線程數(shù)
phred<qual>： 設(shè)定堿基質(zhì)量編碼模式，兩種模式 `-phred33` 或 `-phred64`股缸， 如果未指定質(zhì)量編碼衡楞，則將自動(dòng)識(shí)別, `fastq` 質(zhì)量值編碼可以參考：[https://wiki.bits.vib.be/index.php/FASTQ](https://wiki.bits.vib.be/index.php/FASTQ)
ILUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>： 該步驟用于尋找并去除Illumina接頭
<fastaWithAdaptersEtc1>: 指定包含所有接頭、PCR序列等的fasta文件的路徑;
<seed mismatches>： 指定允許執(zhí)行完全匹配的最大不匹配數(shù)敦姻；
<palindrome clip threshold>：指定兩個(gè)成對(duì)接頭reads之間的匹配對(duì)于雙端回文read比對(duì)的精度瘾境；
<Simple clip threshold>： 指定任何接頭等序列與read之間的匹配精度閾值；
<minAdapterLength>: 設(shè)定接頭的最小長(zhǎng)度镰惦。如果未指定迷守，默認(rèn)為8個(gè)bp;
<keepBothReads>： 在回文模式檢測(cè)到read測(cè)穿并刪除接頭序列后，反向讀取包含與正向讀取相同的序列信息旺入。因此兑凿，默認(rèn)行為是完全刪除反向讀取。通過為該參數(shù)指定true茵瘾，反向讀取也將被保留
SLIDINGWINDOW:<windowSize>:<requiredQuality>: 滑窗修剪, 它從5'端開始掃描礼华，當(dāng)窗口內(nèi)的平均質(zhì)量低于閾值時(shí)，剔除該窗口內(nèi)的所有堿基;
<windowSize>: 設(shè)定窗口大小, 覆蓋堿基數(shù)量龄捡；
<requiredQuality>: 設(shè)定平均質(zhì)量卓嫂。
HEADCROP:<length> ：切除read起始端低于閾值的堿基
<length> 從`read` 起始端開始要切除的長(zhǎng)度
TRAILING:<quality>：從末端移除低質(zhì)量的堿基。只要堿基的質(zhì)量值低于閾值聘殖，則切除該堿基晨雳，并調(diào)查下一個(gè)堿基（因?yàn)門rimmomatic從3'primeend開始，將是位于剛切除堿基之前的堿基）奸腺。 此方法可用于去除Illumina低質(zhì)量段區(qū)域（質(zhì)量分?jǐn)?shù)標(biāo)記為2）,
<quality>: 指定保留堿基所需的最低質(zhì)量
MINLEN:<length>: 設(shè)置保留reads的最小長(zhǎng)度餐禁。
<length>:可被保留的最短 read 長(zhǎng)度。
trimmomatic PE -threads 4 -phred33 \
raw-data/A1_1.fastq \
raw-data/A1_2.fastq \
data_analysis/trimming/reads/A1_1.fastq \
data_analysis/trimming/reads/A1_1.singleton.fastq \
data_analysis/trimming/reads/A1_2.fastq \
data_analysis/trimming/reads/A1_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36
trimmomatic PE -threads 4 -phred33 \
raw-data/A2_1.fastq \
raw-data/A2_2.fastq \
data_analysis/trimming/reads/A2_1.fastq \
data_analysis/trimming/reads/A2_1.singleton.fastq \
data_analysis/trimming/reads/A2_2.fastq \
data_analysis/trimming/reads/A2_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36
trimmomatic PE -threads 4 -phred33 \
raw-data/D1_1.fastq \
raw-data/D1_2.fastq \
data_analysis/trimming/reads/D1_1.fastq \
data_analysis/trimming/reads/D1_1.singleton.fastq \
data_analysis/trimming/reads/D1_2.fastq \
data_analysis/trimming/reads/D1_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36
trimmomatic PE -threads 4 -phred33 \
raw-data/D2_1.fastq \
raw-data/D2_2.fastq \
data_analysis/trimming/reads/D2_1.fastq \
data_analysis/trimming/reads/D2_1.singleton.fastq \
data_analysis/trimming/reads/D2_2.fastq \
data_analysis/trimming/reads/D2_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3
trimmomatic PE -threads 4 -phred33 \
raw-data/H1_1.fastq \
raw-data/H1_2.fastq \
data_analysis/trimming/reads/H1_1.fastq \
data_analysis/trimming/reads/H1_1.singleton.fastq \
data_analysis/trimming/reads/H1_2.fastq \
data_analysis/trimming/reads/H1_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36
trimmomatic PE -threads 4 -phred33 \
raw-data/H2_1.fastq \
raw-data/H2_2.fastq \
data_analysis/trimming/reads/H2_1.fastq \
data_analysis/trimming/reads/H2_1.singleton.fastq \
data_analysis/trimming/reads/H2_2.fastq \
data_analysis/trimming/reads/H2_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36
將雙端數(shù)據(jù)進(jìn)行合并
創(chuàng)建目錄 :
mkdir -p data_analysis/mergepairs/reads
合并雙端序列 :
使用usearch [-fastq_mergepairs] 命令合并A1的雙端序列突照，并輸出結(jié)果到data_analysis/mergepairs/reads/目錄下:
usearch -fastq_mergepairs參數(shù)解析
-fastq_mergepairs 正向FASTQ文件名（輸入）
-reverse 反向FASTQ文件名（輸入）
-fastqout 合并后的FASTQ文件名（輸出）
-fastqout_notmerged_fwd 未合并的正向序列FASTQ文件名（輸出）
-fastqout_notmerged_rev 未合并的反向序列FASTQ文件名（輸出）
-log log文件名（輸出）
-report 總結(jié)報(bào)告文件名（輸出）
-threads 線程數(shù)
-fastq_minmergelen 合并序列的最小長(zhǎng)度
-fastq_maxmergelen 合并序列的最大長(zhǎng)度
-fastq_maxdiffs 最大錯(cuò)配數(shù)帮非，默認(rèn)為5
-fastq_pctid 最小序列質(zhì)量，默認(rèn)為90%
-fastq_minovlen 合并后的序列如果低于指定值則過濾讹蘑，默認(rèn)為16
-fastq_trunctail 在第一個(gè)小于指定Q值處截?cái)嘈蛄心┛J(rèn)為
-fastq_minlen 如果-fastq_trunctail截?cái)嗪蟮男蛄行∮谥付ㄖ担瑒t過濾掉序列
usearch -fastq_mergepairs data_analysis/trimming/reads/A1_1.fastq \
-reverse data_analysis/trimming/reads/A1_2.fastq \
-fastqout data_analysis/mergepairs/reads/A1.fastq \
-fastqout_notmerged_fwd data_analysis/mergepairs/reads/A1.notmerged_fwd.fastq \
-fastqout_notmerged_rev data_analysis/mergepairs/reads/A1.notmerged_rev.fastq \
-log data_analysis/mergepairs/reads/A1.log \
-report data_analysis/mergepairs/reads/A1.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10 \
-fastq_pctid 90 \
-fastq_minovlen 16 \
-fastq_trunctail 2 \
-fastq_minlen 64
usearch -fastq_mergepairs data_analysis/trimming/reads/A2_1.fastq \
-reverse data_analysis/trimming/reads/A2_2.fastq \
-fastqout data_analysis/mergepairs/reads/A2.fastq \
-fastqout_notmerged_fwd data_analysis/mergepairs/reads/A2.notmerged_fwd.fastq \
-fastqout_notmerged_rev data_analysis/mergepairs/reads/A2.notmerged_rev.fastq \
-log data_analysis/mergepairs/reads/A2.log \
-report data_analysis/mergepairs/reads/A2.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10 \
-fastq_pctid 90 \
-fastq_minovlen 16 \
-fastq_trunctail 2 \
-fastq_minlen 64
usearch -fastq_mergepairs data_analysis/trimming/reads/D1_1.fastq \
-reverse data_analysis/trimming/reads/D1_2.fastq \
-fastqout data_analysis/mergepairs/reads/D1.fastq \
-fastqout_notmerged_fwd data_analysis/mergepairs/reads/D1.notmerged_fwd.fastq \
-fastqout_notmerged_rev data_analysis/mergepairs/reads/D1.notmerged_rev.fastq \
-log data_analysis/mergepairs/reads/D1.log \
-report data_analysis/mergepairs/reads/D1.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10 \
-fastq_pctid 90 \
-fastq_minovlen 16 \
-fastq_trunctail 2 \
-fastq_minlen 64
usearch -fastq_mergepairs data_analysis/trimming/reads/D2_1.fastq \
-reverse data_analysis/trimming/reads/D2_2.fastq \
-fastqout data_analysis/mergepairs/reads/D2.fastq \
-fastqout_notmerged_fwd data_analysis/mergepairs/reads/D2.notmerged_fwd.fastq \
-fastqout_notmerged_rev data_analysis/mergepairs/reads/D2.notmerged_rev.fastq \
-log data_analysis/mergepairs/reads/D2.log \
-report data_analysis/mergepairs/reads/D2.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10 \
-fastq_pctid 90 \
-fastq_minovlen 16 \
-fastq_trunctail 2 \
-fastq_minlen 64
usearch -fastq_mergepairs data_analysis/trimming/reads/H1_1.fastq \
-reverse data_analysis/trimming/reads/H1_2.fastq \
-fastqout data_analysis/mergepairs/reads/H1.fastq \
-fastqout_notmerged_fwd data_analysis/mergepairs/reads/H1.notmerged_fwd.fastq \
-fastqout_notmerged_rev data_analysis/mergepairs/reads/H1.notmerged_rev.fastq \
-log data_analysis/mergepairs/reads/H1.log \
-report data_analysis/mergepairs/reads/H1.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10 \
-fastq_pctid 90 \
-fastq_minovlen 16 \
-fastq_trunctail 2 \
-fastq_minlen 64
usearch -fastq_mergepairs data_analysis/trimming/reads/H2_1.fastq \
-reverse data_analysis/trimming/reads/H2_2.fastq \
-fastqout data_analysis/mergepairs/reads/H2.fastq \
-fastqout_notmerged_fwd data_analysis/mergepairs/reads/H2.notmerged_fwd.fastq \
-fastqout_notmerged_rev data_analysis/mergepairs/reads/H2.notmerged_rev.fastq \
-log data_analysis/mergepairs/reads/H2.log \
-report data_analysis/mergepairs/reads/H2.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10 \
-fastq_pctid 90 \
-fastq_minovlen 16 \
-fastq_trunctail 2 \
-fastq_minlen 64
** 統(tǒng)計(jì)序列信息 :**
使用atlas-utils [fqchk] 統(tǒng)計(jì)序列質(zhì)量
cat data_analysis/mergepairs/reads/A1.fastq \
|atlas-utils \
fqchk -q 33 -l A1 \
- >data_analysis/mergepairs/reads/A1.txt
cat data_analysis/mergepairs/reads/A2.fastq \
|atlas-utils \
fqchk -q 33 -l A2 \
- >data_analysis/mergepairs/reads/A2.txt
cat data_analysis/mergepairs/reads/D1.fastq \
|atlas-utils \
fqchk -q 33 -l D1 \
- >data_analysis/mergepairs/reads/D1.txt
cat data_analysis/mergepairs/reads/D2.fastq \
|atlas-utils \
fqchk -q 33 -l D2 \
- >data_analysis/mergepairs/reads/D2.txt
cat data_analysis/mergepairs/reads/H1.fastq \
|atlas-utils \
fqchk -q 33 -l H1 \
- >data_analysis/mergepairs/reads/H1.txt
cat data_analysis/mergepairs/reads/H2.fastq \
|atlas-utils \
fqchk -q 33 -l H2 \
- >data_analysis/mergepairs/reads/H2.txt
統(tǒng)計(jì)信息 :
合并所有mergepairs/reads/目錄下后綴為.txt的文件并傳遞給tsv-utils view 去除中間的注釋行, 使用usearch-utils mergepairs統(tǒng)計(jì)log文件
mkdir -p data_analysis/mergepairs/report
cat data_analysis/mergepairs/reads/*.txt | tsv-utils view - >data_analysis/mergepairs/report/sequencing.stats.txt
比對(duì)引物信息并去除
創(chuàng)建引物信息目錄
mkdir -p data_analysis/primer_strip/reads
usearch -search_pcr2 data_analysis/mergepairs/reads/A1.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/A1.txt \
-fastqout data_analysis/primer_strip/reads/A1.fastq \
-log data_analysis/primer_strip/reads/A1.log
usearch -search_pcr2 data_analysis/mergepairs/reads/A2.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/A2.txt \
-fastqout data_analysis/primer_strip/reads/A2.fastq \
-log data_analysis/primer_strip/reads/A2.log
usearch -search_pcr2 data_analysis/mergepairs/reads/D1.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/D1.txt \
-fastqout data_analysis/primer_strip/reads/D1.fastq \
-log data_analysis/primer_strip/reads/D1.log
usearch -search_pcr2 data_analysis/mergepairs/reads/H1.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/H1.txt \
-fastqout data_analysis/primer_strip/reads/H1.fastq \
-log data_analysis/primer_strip/reads/H1.log
usearch -search_pcr2 data_analysis/mergepairs/reads/D2.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/D2.txt \
-fastqout data_analysis/primer_strip/reads/D2.fastq \
-log data_analysis/primer_strip/reads/D2.log
usearch -search_pcr2 data_analysis/mergepairs/reads/H2.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/H2.txt \
-fastqout data_analysis/primer_strip/reads/H2.fastq \
-log data_analysis/primer_strip/reads/H2.log
導(dǎo)入數(shù)據(jù)
首先需要構(gòu)建數(shù)據(jù)信息表格
mkdir -p data_analysis/dada2/report
mkdir -p data_analysis/dada2/qza/
用于輸入到QIIME 2的所有數(shù)據(jù)都以QIIME 2對(duì)象的形式出現(xiàn)座慰，其中包含有關(guān)數(shù)據(jù)類型和數(shù)據(jù)源的信息陨舱。因此，我們需要做的第一件事是將這些序列數(shù)據(jù)文件導(dǎo)入到QIIME 2對(duì)象中版仔。
開始導(dǎo)入數(shù)據(jù)
qiim-utils manifest mapping_file.txt \
/project/qiime-tk/larry/data_analysis/primer_strip/reads \
> data_analysis/dada2/report/manifest.txt
這個(gè)QIIME 2對(duì)象的語(yǔ)義類型是SequencesWithQuality游盲。 QIIME 對(duì)象SequencesWithQuality是帶有質(zhì)量序列的樣本误墓，。要導(dǎo)入其他格式的序列數(shù)據(jù)益缎，請(qǐng)參閱導(dǎo)入數(shù)據(jù)教程谜慌。
time qiime tools import --type 'SampleData[SequencesWithQuality]' \
--input-path data_analysis/dada2/report/manifest.txt \
--output-path data_analysis/dada2/qza/demux.qza \
--input-format SingleEndFastqManifestPhred33V2
結(jié)果統(tǒng)計(jì)
time qiime demux summarize \
--i-data data_analysis/dada2/qza/demux.qza \
--o-visualization data_analysis/dada2/qza/demux.qzv
2019.7.4.01.jpg918×1877](http://bailab.genetics.ac.cn/markdown/qiime2/fig/2019.7.4.01.jpg "2019.7.4.01.jpg")
圖1. 樣本拆分結(jié)果統(tǒng)計(jì)結(jié)果——樣本數(shù)據(jù)量可視化圖表。
主要分為三部分：上部為摘要莺奔；中部為樣本不同數(shù)據(jù)量分布頻率柱狀圖欣范，可下載PDF，下部為每個(gè)樣本的測(cè)序量令哟。上方面板還可切換至交互式質(zhì)量圖Interactive Qaulity Plot頁(yè)面熙卡。如下圖2。（此處再次點(diǎn)贊??QIIME 2的強(qiáng)大交互性励饵，幾乎所有的可視化文件都可以最快速度最低成本形成發(fā)表級(jí)可視化文件，大大節(jié)約時(shí)間）
2019.7.4.02.gif838×1412](http://bailab.genetics.ac.cn/markdown/qiime2/fig/2019.7.4.02.gif "2019.7.4.02.gif")
圖2. 交互式質(zhì)量圖Interactive Qaulity Plot查看頁(yè)面滑燃。
同樣為三部分：上部為每個(gè)位置堿基的質(zhì)量分布交互式箱線圖役听，鼠標(biāo)懸停在上面，即可在下面(中部)文字和表格中顯示鼠標(biāo)所在位置堿基質(zhì)量的詳細(xì)信息表窘；下部為拆分樣本的長(zhǎng)度摘要(一般等長(zhǎng)測(cè)序無差別)典予。
qiime tools view demux.qzv
直接在服務(wù)器操作的話可以直接利用這條命令可視化，這條命令的顯示需要圖形界面的支持乐严，如在有圖型界面的Linux上瘤袖，但僅使用SSH登陸方式無法顯示圖形。
版本區(qū)別提示：相比QIIME 1.9.1
在QIIME 1.9.1中昂验，我們一般建議通過QIIME執(zhí)行樣本拆分（例如捂敌，使用split_libraries.py或split_libraries_fastq.py這兩個(gè)腳本），因?yàn)檫@個(gè)步驟還執(zhí)行序列的質(zhì)量控制??〖惹伲現(xiàn)在我們將樣本拆分和質(zhì)量控制步驟分開占婉，因此你可以使用混合多樣本序列（如我們?cè)诖怂龅模┗虿鸱趾蟮男蛄虚_始QIIME 2分析。
在樣本拆分之后甫恩，生成拆分結(jié)果的統(tǒng)計(jì)信息非常重要逆济。統(tǒng)計(jì)結(jié)果允許我們確定每個(gè)樣本獲得多少序列，并且還可以獲得序列數(shù)據(jù)中每個(gè)位置處序列質(zhì)量分布的摘要磺箕。
推薦使用 https://view.qiime2.org 網(wǎng)址顯示結(jié)果
可選使用XShell+XManager支持SSH方式的圖型界面奖慌、虛擬機(jī)圖形界面下或服務(wù)器遠(yuǎn)程桌面方式支持上面命令的圖形結(jié)果。我們上學(xué)期的Linux相關(guān)知識(shí)介紹里面有講到過這個(gè)相關(guān)知識(shí)松靡。以后有機(jī)會(huì)的話简僧，會(huì)在宏基因組流程化數(shù)據(jù)分析里面再講一下。
目前命令行方式想要查看結(jié)果可能很多使用服務(wù)器人員無法實(shí)現(xiàn) (即依賴服務(wù)器安裝了桌面击困，本地依賴XShell+XManager或其它ssh終端和圖形界面軟件)
本地查看可解壓.qzv涎劈，目錄中的data目錄包括詳細(xì)的圖表文件广凸，主要關(guān)注 pdf 和 html 文件，目錄結(jié)構(gòu)如下蛛枚。
── demux
└── 8743ab13-72ca-4adf-9b6c-d97e2dbe8ee3
├── checksums.md5
├── data
│ ├── data.jsonp
│ ├── demultiplex-summary.pdf
│ ├── demultiplex-summary.png
│ ├── dist
│ │ ├── bundle.js
│ │ ├── d3-license.txt
│ │ └── vendor.bundle.js
│ ├── forward-seven-number-summaries.csv
│ ├── index.html
│ ├── overview.html
│ ├── per-sample-fastq-counts.csv
│ ├── q2templateassets
│ │ ├── css
│ │ │ ├── bootstrap.min.css
│ │ │ ├── normalize.css
│ │ │ └── tab-parent.css
│ │ ├── fonts
│ │ │ ├── glyphicons-halflings-regular.eot
│ │ │ ├── glyphicons-halflings-regular.svg
│ │ │ ├── glyphicons-halflings-regular.ttf
│ │ │ ├── glyphicons-halflings-regular.woff
│ │ │ └── glyphicons-halflings-regular.woff2
│ │ ├── img
│ │ │ └── qiime2-rect-200.png
│ │ └── js
│ │ ├── bootstrap.min.js
│ │ ├── child.js
│ │ ├── jquery-3.2.0.min.js
│ │ └── parent.js
│ └── quality-plot.html
├── metadata.yaml
├── provenance
│ ├── action
│ │ └── action.yaml
│ ├── artifacts
│ │ ├── 9594ef07-c414-4658-9345-c726de100d8d
│ │ │ ├── action
│ │ │ │ └── action.yaml
│ │ │ ├── citations.bib
│ │ │ ├── metadata.yaml
│ │ │ └── VERSION
│ │ └── a7a882f3-5e4f-4b5e-8a35-6a1098d21608
│ │ ├── action
│ │ │ ├── action.yaml
│ │ │ └── barcodes.tsv
│ │ ├── citations.bib
│ │ ├── metadata.yaml
│ │ └── VERSION
│ ├── citations.bib
│ ├── metadata.yaml
│ └── VERSION
└── VERSION
qzv文件解壓后文件詳細(xì)谅海，可直接訪問data/index.html打開結(jié)果報(bào)告式網(wǎng)頁(yè)。里面的重要結(jié)果蹦浦，全部可以通過此網(wǎng)頁(yè)進(jìn)行索引扭吁。總而言之最重要的文件其實(shí)就是pdf格式文件和html格式文件盲镶，這集中文件是我們?nèi)粘９ぷ髦姓嬲龝?huì)使用到的文件侥袜。
序列質(zhì)控和生成特征表
Sequence quality control and feature table construction
QIIME 2插件多種質(zhì)量控制方法可選，包括DADA2溉贿、Deblur和基于基本質(zhì)量分?jǐn)?shù)的過濾枫吧。在本教程中，我們使用DADA2和Deblur兩種方法分別介紹這個(gè)步驟宇色。這些步驟是可互相替換的九杂，因此你可以使用自己喜歡的方法。這兩種方法的結(jié)果將是一個(gè)QIIME 2特征表FeatureTable[Frequency]和一個(gè)代表性序列FeatureData[Sequence]對(duì)象宣蠕，Frequency對(duì)象包含數(shù)據(jù)集中每個(gè)樣本中每個(gè)唯一序列的計(jì)數(shù)（頻率）例隆，Sequence對(duì)象將FeatureTable中的特征ID與序列對(duì)應(yīng)。
提醒?：此步主要有DADA2和Deblur兩種方法可選抢蚀，推薦使用DADA2镀层，2016年發(fā)表在Nature Method上，在陰道菌群研究中比OTU聚類結(jié)果看到更多細(xì)節(jié)皿曲；相較USEARCH的UPARSE算法唱逢，目前DADA2方法僅去噪去嵌合，不再按相似度聚類屋休，結(jié)果與真實(shí)物種的序列更接近惶我。
注意??：本節(jié)中此次存在兩種可選方法時(shí)，將創(chuàng)建具有特定方法名稱的對(duì)象（例如博投，使用dada2去噪生成的特性表將被稱為
table.qza）绸贡。在創(chuàng)建這些對(duì)象之后，你將把兩個(gè)選項(xiàng)之一的對(duì)象重命名為更通用的文件名（例如毅哗，table.qza）听怕。為對(duì)象創(chuàng)建特定名稱，然后對(duì)其進(jìn)行重命名的過程僅允許你選擇在本步驟中使用的兩個(gè)選項(xiàng)中之一完成案例虑绵，而不必再次關(guān)注該選項(xiàng)尿瞭。需要注意的是，在這個(gè)步驟或QIIME 2中的任何步驟中翅睛，你給對(duì)象或可視化的文件命名并不重要声搁。
方法1. DADA2
Option 1: DADA2
DADA2是用于檢測(cè)和校正（如果有可能的話）Illumina擴(kuò)增序列數(shù)據(jù)的工作流程黑竞。正如在q2-dada2插件中實(shí)現(xiàn)的，這個(gè)質(zhì)量控制過程將過濾掉在測(cè)序數(shù)據(jù)中鑒定的任何phiX序列（通常存在于標(biāo)記基因Illumina測(cè)序數(shù)據(jù)中疏旨，用于提高擴(kuò)增子測(cè)序質(zhì)量）很魂，并同時(shí)過濾嵌合序列。
DADA2是Susan P. Holmes團(tuán)隊(duì)于2016年發(fā)表于Nature Methods的文章檐涝，截止18年12月22號(hào)Google學(xué)術(shù)統(tǒng)計(jì)引用483次遏匆；DADA2自身也是一套在R語(yǔ)言中完整的擴(kuò)增子分析流程
dada2 denoise-single方法需要兩個(gè)用于質(zhì)量過濾的參數(shù)：--p-trim-left m，它去除每個(gè)序列的前m個(gè)堿基(如引物谁榜、標(biāo)簽序列barcode)幅聘；--p-trunc-len n，它在位置n截?cái)嗝總€(gè)序列窃植。這允許用戶去除序列的低質(zhì)量區(qū)域帝蒿、引物或標(biāo)簽序列等。為了確定要為這兩個(gè)參數(shù)傳遞什么值巷怜，你應(yīng)該查看上面由qiime demux summarize生成的demux.qzv文件中的交互質(zhì)量圖選項(xiàng)卡陵叽。
單端序列去噪, 輸入樣本拆分后結(jié)果；去除左端 0 bp (–p-trim-left丛版，有時(shí)用于切除低質(zhì)量序列、barocde或引物)偏序，序列切成 120 bp 長(zhǎng)(–p-trunc-len)页畦；生成代表序列、特征表和去噪過程統(tǒng)計(jì)研儒。
下面的步驟計(jì)算量較大豫缨。
time qiime dada2 denoise-single --i-demultiplexed-seqs data_analysis/dada2/qza/demux.qza \
--p-trim-left 0 \
--p-trunc-len 0 \
--p-n-threads 0 \
--o-representative-sequences data_analysis/dada2/qza/rep-seqs.qza \
--o-table data_analysis/dada2/qza/table.qza \
--o-denoising-stats data_analysis/dada2/qza/stats.qza
# 實(shí)際計(jì)算時(shí)間，即受服務(wù)器配置影響端朵，還受同臺(tái)服務(wù)器上任務(wù)量影響
生成三個(gè)輸出文件：
-
/stats.qza: dada2計(jì)算統(tǒng)計(jì)結(jié)果好芭。 -
table.qza: 特征表。 -
rep-seqs.qza: 代表序列冲呢。
對(duì)特征表統(tǒng)計(jì)進(jìn)行進(jìn)行可視化
qiime metadata tabulate \
--m-input-file data_analysis/dada2/qza/stats.qza \
--o-visualization data_analysis/dada2/qza/stats.qzv
qiime tools export --input-path data_analysis/dada2/qza/table.qza \
--output-path data_analysis/dada2/denoise ;
qiime tools export --input-path data_analysis/dada2/qza/rep-seqs.qza \
--output-path data_analysis/dada2/denoise ;
qiime tools export --input-path data_analysis/dada2/qza/stats.qza \
--output-path data_analysis/dada2/denoise ;
biom convert \
--input-fp data_analysis/dada2/denoise/feature-table.biom \
--output-fp data_analysis/dada2/denoise/biom-table.tsv --to-tsv;
tail -n +2 data_analysis/dada2/denoise/biom-table.tsv \
|tsv-utils view -r \
- > data_analysis/dada2/denoise/feature-table.tsv
方法2. Deblur（個(gè)人不是很推薦舍败，內(nèi)存占用較大，時(shí)間慢敬拓，所以不演示了邻薯，看一下流程吧）
Deblur使用序列錯(cuò)誤配置文件將錯(cuò)誤的序列與從其來源的真實(shí)生物序列相關(guān)聯(lián)，從而得到高質(zhì)量的序列變異數(shù)據(jù)乘凸，主要為兩個(gè)步驟厕诡。首先，應(yīng)用基于質(zhì)量分?jǐn)?shù)的初始質(zhì)量過濾過程营勤，是Bokulich等人2013年發(fā)表的質(zhì)量過濾方法灵嫌。
Deblur是本軟件作者作為通訊作者2013發(fā)表于Nature Methods的重要擴(kuò)增子代表序列鑒定方法壹罚，截止19年8月25號(hào)Google學(xué)術(shù)統(tǒng)計(jì)引用1259次，
按測(cè)序堿基質(zhì)量過濾序列
time qiime quality-filter q-score \
--i-demux demux.qza \
--o-filtered-sequences demux-filtered.qza \
--o-filter-stats demux-filter-stats.qza
# 用時(shí)40s
輸出對(duì)象:
-
demux-filtered.qza: 序列質(zhì)量過濾后結(jié)果寿羞。 查看 | 下載 -
demux-filter-stats.qza: 序列質(zhì)量過濾后結(jié)果統(tǒng)計(jì)猖凛。 查看 | 下載
接下來，使qiime deblur denoise-16S方法應(yīng)用于Deblur工作流程稠曼。此方法需要一個(gè)用于質(zhì)量過濾的參數(shù)形病，即截?cái)辔恢胣長(zhǎng)度的序列的--p-trim-length n。通常霞幅，Deblur開發(fā)人員建議將該值設(shè)置為質(zhì)量分?jǐn)?shù)中位數(shù)開始下降至低質(zhì)量區(qū)時(shí)的長(zhǎng)度漠吻。在本次數(shù)據(jù)上，質(zhì)量圖（在質(zhì)量過濾之前）表明合理的選擇是在115至130序列位置范圍內(nèi)司恳。這是一個(gè)主觀的評(píng)估途乃。你可能不采用該建議的一種原因是存在多個(gè)批次測(cè)序的元分析。在這種情況的元分析中扔傅，比較所有批次的序列長(zhǎng)度是否相同耍共，以避免人為引入特定的偏差，全局考慮這些是非常重要的猎塞。由于我們已經(jīng)使用修剪長(zhǎng)度為120 bp用于qiime dada2 denoise-single分析试读，并且由于120 bp是基于質(zhì)量圖的結(jié)果，這里我們將使用--p-trim-length 120參數(shù)荠耽。下一個(gè)命令可能需要10分鐘才能運(yùn)行完成钩骇。
deblur去噪16S過程，輸入文件為質(zhì)控后的序列铝量，設(shè)置截取長(zhǎng)度參數(shù)倘屹，生成結(jié)果文件有代表序列、特征表慢叨、樣本統(tǒng)計(jì)纽匙。
time qiime deblur denoise-16S \
--i-demultiplexed-seqs demux-filtered.qza \
--p-trim-length 120 \
--o-representative-sequences rep-seqs-deblur.qza \
--o-table table-deblur.qza \
--p-sample-stats \
--o-stats deblur-stats.qza
-
deblur-stats.qza:過程統(tǒng)計(jì)。 查看 | 下載 -
table-deblur.qza: 特征表拍谐。 查看 | 下載 -
rep-seqs-deblur.qza: 代表序列烛缔。 查看 1 | 下載
time qiime metadata tabulate \
--m-input-file demux-filter-stats.qza \
--o-visualization demux-filter-stats.qzv
time qiime deblur visualize-stats \
--i-deblur-stats deblur-stats.qza \
--o-visualization deblur-stats.qzv
輸出結(jié)果:
示例如下：包括6列力穗，第一列為樣本名稱，2-6列分別為總輸入讀長(zhǎng)气嫁、總保留高讀長(zhǎng)当窗、截?cái)嗟淖x長(zhǎng)、截?cái)嗪筇痰淖x長(zhǎng)和超過最大模糊堿基的讀長(zhǎng)的數(shù)量統(tǒng)計(jì)寸宵。我們通常只關(guān)注2崖面，3列數(shù)量即可元咙，其它列常用于異常的輸助判斷。
| sample-id | total-input-reads | total-retained-reads | reads-truncated | reads-too-short-after-truncation | reads-exceeding-maximum-ambiguous-bases |
|---|---|---|---|---|---|
| #q2:types | numeric | numeric | numeric | numeric | numeric |
| L1S105 | 11340 | 9232 | 10782 | 2066 | 42 |
| L1S140 | 9738 | 8585 | 9459 | 1113 | 40 |
| L1S208 | 11337 | 10149 | 10668 | 1161 | 27 |
2019.7.4.03.jpg1870×760](http://bailab.genetics.ac.cn/markdown/qiime2/fig/2019.7.4.03.jpg "2019.7.4.03.jpg")
圖3. deblur去噪和鑒定ASV處理過程統(tǒng)計(jì)結(jié)果
如果你想用此處結(jié)果下游分析庶香，可以改名為下游分析的起始名稱：
這處演示不運(yùn)行下面兩行代碼，前面添加"#"號(hào)代表注釋简识，需要運(yùn)行請(qǐng)自行刪除行首的“#”
#mv rep-seqs-deblur.qza rep-seqs.qza
#mv table-deblur.qza table.qza
??統(tǒng)一完整的數(shù)據(jù)命名也是非常好的服務(wù)器運(yùn)行習(xí)慣赶掖。
物種組成分析
Taxonomic analysis
在這一節(jié)中，我們將開始探索樣本的物種組成七扰，并將其與樣本元數(shù)據(jù)再次組合奢赂。這個(gè)過程的第一步是為FeatureData[Sequence]的序列進(jìn)行物種注釋。我們將使用經(jīng)過Naive Bayes分類器預(yù)訓(xùn)練的颈走，并由q2-feature-classifier插件來完成這項(xiàng)工作膳灶。這個(gè)分類器是在silva-138-99上訓(xùn)練的，其中序列被修剪到僅包括來自16S區(qū)域的250個(gè)堿基立由，該16S區(qū)域在該分析中采用V4區(qū)域的515F/806R引物擴(kuò)增并測(cè)序轧钓。我們將把這個(gè)分類器應(yīng)用到序列中，并且可以生成從序列到物種注釋結(jié)果關(guān)聯(lián)的可視化锐膜。
注意：物種分類器根據(jù)你特定的樣品制備和測(cè)序參數(shù)進(jìn)行訓(xùn)練時(shí)表現(xiàn)最好毕箍，包括用于擴(kuò)增的引物和測(cè)序序列的長(zhǎng)度。因此道盏，一般來說而柑，你應(yīng)該按照使用
q2-feature-classifier的訓(xùn)練特征分類器的說明來訓(xùn)練自己的物種分類器。我們?cè)?a target="_blank">數(shù)據(jù)資源頁(yè)面上提供了一些通用的分類器捞奕，包括基于Silva的16S分類器，不過將來我們可能會(huì)停止提供這些分類器拄轻，而讓用戶訓(xùn)練他們自己的分類器颅围，這將與他們的序列數(shù)據(jù)最相關(guān)。
mkdir -p data_analysis/classify/sklearn
提醒?：這一步會(huì)占用大量線程所以請(qǐng)大家不要操作恨搓，我這邊演示后會(huì)提供在上層目錄中
qiime feature-classifier classify-sklearn --i-classifier ../db/silva-138-99-nb-classifier.qza \
--i-reads data_analysis/dada2/qza/rep-seqs.qza \
--p-n-jobs 4 \
--o-classification data_analysis/classify/sklearn/taxonomy.qza
qiime metadata tabulate \
--m-input-file data_analysis/classify/sklearn/taxonomy.qza \
--o-visualization data_analysis/classify/sklearn/taxonomy.qzv
我們可以在網(wǎng)頁(yè)端口導(dǎo)入可視化文件來看一下具體的結(jié)果
導(dǎo)出物種注釋信息
qiime tools export --input-path data_analysis/classify/sklearn/taxonomy.qzv \
--output-path data_analysis/classify/taxonomy
基于此我們就得到了樣品中的微生物信息可以了解樣品中究竟有哪些物種在其中
qiime taxa barplot \
--i-table data_analysis/dada2/qza/table.qza \
--i-taxonomy data_analysis/classify/sklearn/taxonomy.qza \
--m-metadata-file mapping_file.txt \
--o-visualization data_analysis/classify/sklearn/taxa-bar-plots.qzv
qiime tools export \
--input-path data_analysis/classify/sklearn/taxa-bar-plots.qzv \
--output-path data_analysis/classify/barplot
構(gòu)建進(jìn)化樹用于多樣性分析
Generate a tree for phylogenetic diversity analyses
QIIME 2支持幾種系統(tǒng)發(fā)育多樣性度量方法院促，包括Faith’s Phylogenetic Diversity、weighted和unweighted UniFrac斧抱。除了每個(gè)樣本的特征計(jì)數(shù)（即QIIME2對(duì)象FeatureTable[Frequency]）之外常拓，這些度量還需要將特征彼此關(guān)聯(lián)結(jié)合有根進(jìn)化樹。此信息將存儲(chǔ)在一個(gè)QIIME 2對(duì)象的有根系統(tǒng)發(fā)育對(duì)象Phylogeny[Rooted]中辉浦。為了生成系統(tǒng)發(fā)育樹??弄抬，我們將使用q2-phylogeny插件中的align-to-tree-mafft-fasttree工作流程。
首先宪郊，工作流程使用mafft程序執(zhí)行對(duì)FeatureData[Sequence]中的序列進(jìn)行多序列比對(duì)掂恕，以創(chuàng)建QIIME 2對(duì)象FeatureData[AlignedSequence]拖陆。接下來，流程屏蔽（mask或過濾）對(duì)齊的的高度可變區(qū)(高變區(qū))懊亡，這些位置通常被認(rèn)為會(huì)增加系統(tǒng)發(fā)育樹的噪聲依啰。隨后，流程應(yīng)用FastTree基于過濾后的比對(duì)結(jié)果生成系統(tǒng)發(fā)育樹店枣。FastTree程序創(chuàng)建的是一個(gè)無根樹速警，因此在本節(jié)的最后一步中，應(yīng)用根中點(diǎn)法將樹的根放置在無根樹中最長(zhǎng)端到端距離的中點(diǎn)鸯两，從而形成有根樹闷旧。
mkdir -p data_analysis/phylogeny/qza
qiime phylogeny align-to-tree-mafft-fasttree --i-sequences data_analysis/dada2/qza/rep-seqs.qza \
--o-alignment data_analysis/phylogeny/qza/aligned-rep-seqs.qza \
--o-masked-alignment data_analysis/phylogeny/qza/masked-aligned-rep-seqs.qza \
--o-tree data_analysis/phylogeny/qza/unrooted-tree.qza \
--o-rooted-tree data_analysis/phylogeny/qza/rooted-tree.qza
qiime tools export --input-path data_analysis/phylogeny/qza/rooted-tree.qza \
--output-path data_analysis/phylogeny/tree
多樣性分析
根據(jù)最小測(cè)序深度計(jì)算物種多樣性信息
rm -rf data_analysis/metrics/diversity
qiime diversity core-metrics-phylogenetic \
--i-phylogeny data_analysis/phylogeny/qza/rooted-tree.qza \
--i-table data_analysis/dada2/qza/table.qza \
--p-sampling-depth 6000 \
--m-metadata-file mapping_file.txt \
--output-dir data_analysis/metrics/diversity
獲得特征序列表
qiime tools export --input-path data_analysis/metrics/diversity/rarefied_table.qza \
--output-path data_analysis/metrics/report ;
biom convert --input-fp data_analysis/metrics/report/feature-table.biom \
--output-fp data_analysis/metrics/report/biom-table.tsv --to-tsv;
tail -n +2 data_analysis/metrics/report/biom-table.tsv \
|tsv-utils view -r \
- >data_analysis/metrics/report/feature_table.tsv
Alpha和beta多樣性分析
Alpha and beta diversity analysis
此處多樣性分析的基礎(chǔ)知識(shí)和相關(guān)方法眾多，詳細(xì)參考我在公眾號(hào)平臺(tái)上的總結(jié)甩卓，此處不論述鸠匀。QIIME 2的多樣性分析使用q2-diversity插件，該插件支持計(jì)算α和β多樣性指數(shù)逾柿、并應(yīng)用相關(guān)的統(tǒng)計(jì)檢驗(yàn)以及生成交互式可視化圖表缀棍。我們將首先應(yīng)用core-metrics-phylogenetic方法，該方法將FeatureTable[Frequency](特征表[頻率])抽平到用戶指定的測(cè)序深度机错，然后計(jì)算幾種常用的α和β多樣性指數(shù)爬范，并使用Emperor為每個(gè)β多樣性指數(shù)生成主坐標(biāo)分析（PCoA）圖。默認(rèn)情況下計(jì)算的方法有：
劃重點(diǎn)：理解下面4種alpha和beta多樣性指數(shù)的所代表的生物學(xué)意義至關(guān)重要??弱匪。
- α多樣性
- 香農(nóng)(Shannon’s)多樣性指數(shù)（群落豐富度的定量度量青瀑，即包括豐富度
richness和均勻度evenness兩個(gè)層面） - 可觀測(cè)的OTU(Observed OTUs，群落豐富度的定性度量萧诫，只包括豐富度）
- Faith’s系統(tǒng)發(fā)育多樣性（包含特征之間的系統(tǒng)發(fā)育關(guān)系的群落豐富度的定性度量）
- 均勻度Evenness（或 Pielou’s均勻度像鸡；群落均勻度的度量）
- 香農(nóng)(Shannon’s)多樣性指數(shù)（群落豐富度的定量度量青瀑，即包括豐富度
- β多樣性
- Jaccard距離（群落差異的定性度量，即只考慮種類周崭，不考慮豐度）
- Bray-Curtis距離（群落差異的定量度量撒璧，較常用）
- 非加權(quán)UniFrac距離（包含特征之間的系統(tǒng)發(fā)育關(guān)系的群落差異定性度量）
- 加權(quán)UniFrac距離（包含特征之間的系統(tǒng)發(fā)育關(guān)系的群落差異定量度量）
需要提供給這個(gè)腳本的一個(gè)重要參數(shù)是--p-sampling-depth，它是指定重采樣（即稀疏/稀疏rarefaction）深度及刻。因?yàn)榇蠖鄶?shù)多樣指數(shù)對(duì)不同樣本的不同測(cè)序深度敏感镀裤，所以這個(gè)腳本將隨機(jī)地將每個(gè)樣本的測(cè)序量重新采樣至該參數(shù)值。例如缴饭，提供--p-sampling-depth 500暑劝，則此步驟將對(duì)每個(gè)樣本中的計(jì)數(shù)進(jìn)行無放回抽樣，從而使得結(jié)果表中的每個(gè)樣本的總計(jì)數(shù)為500颗搂。如果任何樣本的總計(jì)數(shù)小于該值担猛，那么這些樣本將從多樣性分析中刪除。選擇這個(gè)值很棘手。建議通過查看上面創(chuàng)建的表table.qzv文件中呈現(xiàn)的信息并選擇一個(gè)盡可能高的值（因此每個(gè)樣本保留更多的序列）同時(shí)盡可能少地排除樣本來進(jìn)行選擇毁习。
下面多樣性分析智嚷，需要基于重采樣/抽平(rarefaction)標(biāo)準(zhǔn)化的特征表，標(biāo)準(zhǔn)化采用無放回重抽樣至序列一致纺且，如何設(shè)計(jì)樣品重采樣深度參數(shù)
--p-sampling-depth呢盏道？
如是數(shù)據(jù)量都很大，選最小的即可载碌。如果有個(gè)別數(shù)據(jù)量非常小猜嘱，去除最小值再選最小值。比如此分析最小值為917嫁艇，我們選擇1109深度重采樣朗伶，即保留了大部分樣品用于分析，又去除了數(shù)據(jù)量過低的異常值步咪。本示例為近10年前測(cè)序技術(shù)的通量水平论皆，454測(cè)序時(shí)代抽平至1000條即可，現(xiàn)在看來數(shù)據(jù)量很小猾漫。目錄一般采用HiSeq2500或NovaSeq6000的 PE250模式測(cè)序点晴，數(shù)據(jù)量都非常大，通趁踔埽可以采用3萬或5萬的標(biāo)準(zhǔn)抽平粒督，仍可保留90%以上樣本。過低或過高一般結(jié)果也會(huì)波動(dòng)較大禽翼，不建議放在一起分析屠橄。
qiime tools export \
--input-path data_analysis/metrics/diversity/shannon_vector.qza \
--output-path data_analysis/alpha/diversity ;
mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/shannon.tsv ;
qiime tools export \
--input-path data_analysis/metrics/diversity/faith_pd_vector.qza \
--output-path data_analysis/alpha/diversity ;
mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/faith_pd.tsv ;
qiime tools export \
--input-path data_analysis/metrics/diversity/observed_features_vector.qza \
--output-path data_analysis/alpha/diversity ;
mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/observed_features.tsv ;
qiime tools export \
--input-path data_analysis/metrics/diversity/evenness_vector.qza \
--output-path data_analysis/alpha/diversity ;
mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/evenness.tsv
計(jì)算稀釋度曲線
mkdir -p data_analysis/alpha/qza
qiime diversity alpha-rarefaction \
--p-min-depth 1 \
--p-steps 10 \
--p-iterations 10 \
--p-metrics shannon simpson chao1 observed_features goods_coverage faith_pd \
--i-table data_analysis/dada2/qza/table.qza \
--i-phylogeny data_analysis/phylogeny/qza/rooted-tree.qza \
--p-max-depth 7500 \
--m-metadata-file mapping_file.txt \
--o-visualization data_analysis/alpha/qza/alpha-rarefaction.qzv ;
qiime tools export \
--input-path data_analysis/alpha/qza/alpha-rarefaction.qzv \
--output-path data_analysis/alpha/rarefaction ;
beta多樣性
mkdir -p data_analysis/beta/matrix
qiime tools export \
--input-path data_analysis/metrics/diversity/jaccard_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/jaccard.tsv ;
qiime tools export --input-path data_analysis/metrics/diversity/bray_curtis_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/bray_curtis.tsv ;
qiime tools export \
--input-path data_analysis/metrics/diversity/weighted_unifrac_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/weighted_unifrac.tsv ;
qiime tools export \
--input-path data_analysis/metrics/diversity/unweighted_unifrac_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/unweighted_unifrac.tsv
mkdir -p data_analysis/beta/emperor/jaccard ;
qiime tools export \
--input-path data_analysis/metrics/diversity/jaccard_emperor.qzv \
--output-path data_analysis/beta/emperor/jaccard ;
mkdir -p data_analysis/beta/emperor/bray_curtis ;
qiime tools export \
--input-path data_analysis/metrics/diversity/bray_curtis_emperor.qzv \
--output-path data_analysis/beta/emperor/bray_curtis ;
mkdir -p data_analysis/beta/emperor/weighted_unifrac ;
qiime tools export \
--input-path data_analysis/metrics/diversity/weighted_unifrac_emperor.qzv \
--output-path data_analysis/beta/emperor/weighted_unifrac ;
mkdir -p data_analysis/beta/emperor/unweighted_unifrac ;
qiime tools export --input-path data_analysis/metrics/diversity/unweighted_unifrac_emperor.qzv \
--output-path data_analysis/beta/emperor/unweighted_unifrac
PCoA多樣性分析
mkdir -p data_analysis/beta/pcoa/jaccard ;
qiime tools export \
--input-path data_analysis/metrics/diversity/jaccard_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/jaccard ;
mkdir -p data_analysis/beta/pcoa/bray_curtis ;
qiime tools export \
--input-path data_analysis/metrics/diversity/bray_curtis_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/bray_curtis ;
mkdir -p data_analysis/beta/pcoa/weighted_unifrac ;
qiime tools export \
--input-path data_analysis/metrics/diversity/weighted_unifrac_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/weighted_unifrac ;
mkdir -p data_analysis/beta/pcoa/unweighted_unifrac ;
qiime tools export \
--input-path data_analysis/metrics/diversity/unweighted_unifrac_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/unweighted_unifrac
rank_abundance
mkdir -p data_analysis/rank_abundance/report/
rank_abundance.R data_analysis/metrics/report/feature_table.tsv \
data_analysis/rank_abundance/report/rank_abundance.pdf ;
pdf2png data_analysis/rank_abundance/report/rank_abundance.pdf ;
specaccum_curve
mkdir -p data_analysis/specaccum_curve/report/
specaccum_curve.R data_analysis/metrics/report/feature_table.tsv \
data_analysis/specaccum_curve/report/specaccum_curve.pdf
pdf2png data_analysis/specaccum_curve/report/specaccum_curve.pdf
picrust2 功能預(yù)測(cè)
這一步也需要占據(jù)大量的計(jì)算資源所以請(qǐng)大家不要自行分析，請(qǐng)看我進(jìn)行演示
deactive
conda activate /biostack/bioconda/envs/picrust2
rm -rf data_analysis/picrust2/pipeline;
picrust2_pipeline.py \
-s data_analysis/dada2/denoise/dna-sequences.fasta \
-i data_analysis/metrics/report/feature_table.tsv \
-o data_analysis/picrust2/pipeline \
-p 64 --no_pathways \
--stratified \
--wide_table \
--per_sequence_contrib
mkdir -p data_analysis/picrust2/prediction
mkdir -p data_analysis/picrust2/report
tsv-utils strip mapping_file.txt | cut -f1,2 | grep -v '#' >data_analysis/picrust2/report/metadata.txt;
tsv-utils definition \
-d " " /biostack/tools/protocols/qiim-tk-0.0.4/bins/../db/picrust2/ko_info.tsv data_analysis/picrust2/pipeline/KO_metagenome_out/pred_metagenome_unstrat.tsv.gz > data_analysis/picrust2/prediction/ko.txt ;
tsv-utils definition \
-d " " /biostack/tools/protocols/qiim-tk-0.0.4/bins/../db/picrust2/ec_level4_info.tsv data_analysis/picrust2/pipeline/EC_metagenome_out/pred_metagenome_unstrat.tsv.gz \
>data_analysis/picrust2/prediction/enzyme.txt
功能預(yù)測(cè)基礎(chǔ)可視化
qiime2
mkdir -p data_analysis/picrust2/prediction/ko ;
mkdir -p data_analysis/picrust2/prediction/enzyme
PCA.R -g T data_analysis/picrust2/prediction/ko.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/ko/ko.pca.pdf;
pdf2png data_analysis/picrust2/prediction/ko/ko.pca.pdf ;
PCoA.R -g T -m bray data_analysis/picrust2/prediction/ko.txt data_analysis/picrust2/report/metadata.txt /project/qiime-tk/data_analysis/picrust2/prediction/ko/ko.pcoa.pdf;
pdf2png data_analysis/picrust2/prediction/ko/ko.pcoa.pdf ;
NMDS.R -g T -m bray data_analysis/picrust2/prediction/ko.txt data_analysis/picrust2/report/metadata.txt /project/qiime-tk/data_analysis/picrust2/prediction/ko/ko.nmds.pdf;
pdf2png data_analysis/picrust2/prediction/ko/ko.nmds.pdf ;
mkdir -p data_analysis/picrust2/prediction/enzyme ;
PCA.R -g T data_analysis/picrust2/prediction/enzyme.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/enzyme/enzyme.pca.pdf;
pdf2png data_analysis/picrust2/prediction/enzyme/enzyme.pca.pdf ;
PCoA.R -g T -m bray data_analysis/picrust2/prediction/enzyme.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/enzyme/enzyme.pcoa.pdf;
pdf2png data_analysis/picrust2/prediction/enzyme/enzyme.pcoa.pdf ;
NMDS.R -g T -m bray data_analysis/picrust2/prediction/enzyme.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/enzyme/enzyme.nmds.pdf;
pdf2png data_analysis/picrust2/prediction/enzyme/enzyme.nmds.pdf
######################################################
