內(nèi)部培訓(xùn)資料（2022.10.29）
QIIME 2分析實(shí)例--人類(lèi)腸癌不同部位微生物多樣性分析
“Moving Pictures” tutorial
熊金波實(shí)驗(yàn)室出品 ??????內(nèi)部交流培訓(xùn)使用
寧波大學(xué)海洋學(xué)院 Larry 陸
本節(jié)1.6萬(wàn)字，14張圖。閱讀時(shí)間大約50分鐘。
在本次分享中，我們將使用QIIME2軟件進(jìn)行微生物數(shù)據(jù)分析。在五個(gè)時(shí)間點(diǎn)對(duì)來(lái)自?xún)蓚€(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一直以來(lái)被詬病最多的就是非常難以安裝，目前據(jù)我了解有三種主流安裝方法，根據(jù)環(huán)境任選其一即可，實(shí)際的安裝涉及到服務(wù)器運(yùn)維工作，與我們實(shí)驗(yàn)室日常工作相關(guān)性不大，所以此處不做過(guò)多介紹。我們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ù)分析問(wèn)題。
我們?cè)诿看畏治鲩_(kāi)始前，必須先進(jìn)入工作目錄，除非你是一個(gè)把什么東西都放在桌面上還很工作更有效率的人??。
# 定義工作目錄變量，方便以后多次使用
wd=~/github/QIIME2ChineseManual/2019.7
mkdir -p $wd
# 進(jìn)入工作目錄，是不是很簡(jiǎn)介，這樣無(wú)論你在什么位置就可以快速回到項(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版本較老用戶(hù)，使用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
正式開(kāi)始分析數(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)行為是完全刪除反向讀取。通過(guò)為該參數(shù)指定true，反向讀取也將被保留
SLIDINGWINDOW:<windowSize>:<requiredQuality>: 滑窗修剪, 它從5'端開(kāi)始掃描，當(dāng)窗口內(nèi)的平均質(zhì)量低于閾值時(shí)，剔除該窗口內(nèi)的所有堿基;
<windowSize>: 設(shè)定窗口大小, 覆蓋堿基數(shù)量；
<requiredQuality>: 設(shè)定平均質(zhì)量。
HEADCROP:<length> ：切除read起始端低于閾值的堿基
<length> 從`read` 起始端開(kāi)始要切除的長(zhǎng)度
TRAILING:<quality>：從末端移除低質(zhì)量的堿基。只要堿基的質(zhì)量值低于閾值，則切除該堿基，并調(diào)查下一個(gè)堿基（因?yàn)門(mén)rimmomatic從3'primeend開(kāi)始，將是位于剛切除堿基之前的堿基）。 此方法可用于去除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 合并后的序列如果低于指定值則過(guò)濾，默認(rèn)為16
-fastq_trunctail 在第一個(gè)小于指定Q值處截?cái)嘈蛄校J(rèn)為
-fastq_minlen 如果-fastq_trunctail截?cái)嗪蟮男蛄行∮谥付ㄖ担瑒t過(guò)濾掉序列
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ù)類(lèi)型和數(shù)據(jù)源的信息。因此，我們需要做的第一件事是將這些序列數(shù)據(jù)文件導(dǎo)入到QIIME 2對(duì)象中。
開(kāi)始導(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ǔ)義類(lèi)型是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è)序無(wú)差別)。
qiime tools view demux.qzv
直接在服務(wù)器操作的話可以直接利用這條命令可視化，這條命令的顯示需要圖形界面的支持，如在有圖型界面的Linux上，但僅使用SSH登陸方式無(wú)法顯示圖形。
版本區(qū)別提示：相比QIIME 1.9.1
在QIIME 1.9.1中，我們一般建議通過(guò)QIIME執(zhí)行樣本拆分（例如，使用split_libraries.py或split_libraries_fastq.py這兩個(gè)腳本），因?yàn)檫@個(gè)步驟還執(zhí)行序列的質(zhì)量控制??。現(xiàn)在我們將樣本拆分和質(zhì)量控制步驟分開(kāi)，因此你可以使用混合多樣本序列（如我們?cè)诖怂龅模┗虿鸱趾蟮男蛄虚_(kāi)始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í)介紹里面有講到過(guò)這個(gè)相關(guān)知識(shí)。以后有機(jī)會(huì)的話，會(huì)在宏基因組流程化數(shù)據(jù)分析里面再講一下。
目前命令行方式想要查看結(jié)果可能很多使用服務(wù)器人員無(wú)法實(shí)現(xiàn) (即依賴(lài)服務(wù)器安裝了桌面，本地依賴(lài)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ì)，可直接訪問(wèn)data/index.html打開(kāi)結(jié)果報(bào)告式網(wǎng)頁(yè)。里面的重要結(jié)果，全部可以通過(guò)此網(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ù)的過(guò)濾。在本教程中，我們使用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聚類(lèi)結(jié)果看到更多細(xì)節(jié)；相較USEARCH的UPARSE算法，目前DADA2方法僅去噪去嵌合，不再按相似度聚類(lèi)，結(jié)果與真實(shí)物種的序列更接近。
注意??：本節(jié)中此次存在兩種可選方法時(shí)，將創(chuàng)建具有特定方法名稱(chēng)的對(duì)象（例如，使用dada2去噪生成的特性表將被稱(chēng)為
table.qza）。在創(chuàng)建這些對(duì)象之后，你將把兩個(gè)選項(xiàng)之一的對(duì)象重命名為更通用的文件名（例如，table.qza）。為對(duì)象創(chuàng)建特定名稱(chēng)，然后對(duì)其進(jìn)行重命名的過(guò)程僅允許你選擇在本步驟中使用的兩個(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ì)量控制過(guò)程將過(guò)濾掉在測(cè)序數(shù)據(jù)中鑒定的任何phiX序列（通常存在于標(biāo)記基因Illumina測(cè)序數(shù)據(jù)中，用于提高擴(kuò)增子測(cè)序質(zhì)量），并同時(shí)過(guò)濾嵌合序列。
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ì)量過(guò)濾的參數(shù)：--p-trim-left m，它去除每個(gè)序列的前m個(gè)堿基(如引物、標(biāo)簽序列barcode)；--p-trunc-len n，它在位置n截?cái)嗝總€(gè)序列。這允許用戶(hù)去除序列的低質(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)；生成代表序列、特征表和去噪過(guò)程統(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ò)誤的序列與從其來(lái)源的真實(shí)生物序列相關(guān)聯(lián)，從而得到高質(zhì)量的序列變異數(shù)據(jù)，主要為兩個(gè)步驟。首先，應(yīng)用基于質(zhì)量分?jǐn)?shù)的初始質(zhì)量過(guò)濾過(guò)程，是Bokulich等人2013年發(fā)表的質(zhì)量過(guò)濾方法。
Deblur是本軟件作者作為通訊作者2013發(fā)表于Nature Methods的重要擴(kuò)增子代表序列鑒定方法，截止19年8月25號(hào)Google學(xué)術(shù)統(tǒng)計(jì)引用1259次，
按測(cè)序堿基質(zhì)量過(guò)濾序列
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ì)量過(guò)濾后結(jié)果。 查看 | 下載 -
demux-filter-stats.qza: 序列質(zhì)量過(guò)濾后結(jié)果統(tǒng)計(jì)。 查看 | 下載
接下來(lái)，使qiime deblur denoise-16S方法應(yīng)用于Deblur工作流程。此方法需要一個(gè)用于質(zhì)量過(guò)濾的參數(shù)，即截?cái)辔恢胣長(zhǎng)度的序列的--p-trim-length n。通常，Deblur開(kāi)發(fā)人員建議將該值設(shè)置為質(zhì)量分?jǐn)?shù)中位數(shù)開(kāi)始下降至低質(zhì)量區(qū)時(shí)的長(zhǎng)度。在本次數(shù)據(jù)上，質(zhì)量圖（在質(zhì)量過(guò)濾之前）表明合理的選擇是在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過(guò)程，輸入文件為質(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:過(guò)程統(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é)果:
-
demux-filter-stats.qzv: 質(zhì)量過(guò)程統(tǒng)計(jì)表，同上面提到的stats-dada2.qzv統(tǒng)計(jì)表類(lèi)似。 查看| 下載 1
示例如下：包括6列，第一列為樣本名稱(chēng)，2-6列分別為總輸入讀長(zhǎng)、總保留高讀長(zhǎng)、截?cái)嗟淖x長(zhǎng)、截?cái)嗪筇痰淖x長(zhǎng)和超過(guò)最大模糊堿基的讀長(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處理過(guò)程統(tǒng)計(jì)結(jié)果
如果你想用此處結(jié)果下游分析，可以改名為下游分析的起始名稱(chēng)：
這處演示不運(yùn)行下面兩行代碼，前面添加"#"號(hào)代表注釋?zhuān)枰\(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é)中，我們將開(kāi)始探索樣本的物種組成，并將其與樣本元數(shù)據(jù)再次組合。這個(gè)過(guò)程的第一步是為FeatureData[Sequence]的序列進(jìn)行物種注釋。我們將使用經(jīng)過(guò)Naive Bayes分類(lèi)器預(yù)訓(xùn)練的，并由q2-feature-classifier插件來(lái)完成這項(xiàng)工作。這個(gè)分類(lèi)器是在silva-138-99上訓(xùn)練的，其中序列被修剪到僅包括來(lái)自16S區(qū)域的250個(gè)堿基，該16S區(qū)域在該分析中采用V4區(qū)域的515F/806R引物擴(kuò)增并測(cè)序。我們將把這個(gè)分類(lèi)器應(yīng)用到序列中，并且可以生成從序列到物種注釋結(jié)果關(guān)聯(lián)的可視化。
注意：物種分類(lèi)器根據(jù)你特定的樣品制備和測(cè)序參數(shù)進(jìn)行訓(xùn)練時(shí)表現(xiàn)最好，包括用于擴(kuò)增的引物和測(cè)序序列的長(zhǎng)度。因此，一般來(lái)說(shuō)，你應(yīng)該按照使用
q2-feature-classifier的訓(xùn)練特征分類(lèi)器的說(shuō)明來(lái)訓(xùn)練自己的物種分類(lèi)器。我們?cè)?a target="_blank">數(shù)據(jù)資源頁(yè)面上提供了一些通用的分類(lèi)器，包括基于Silva的16S分類(lèi)器，不過(guò)將來(lái)我們可能會(huì)停止提供這些分類(lèi)器，而讓用戶(hù)訓(xùn)練他們自己的分類(lèi)器，這將與他們的序列數(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)入可視化文件來(lái)看一下具體的結(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)化樹(shù)用于多樣性分析
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)化樹(shù)。此信息將存儲(chǔ)在一個(gè)QIIME 2對(duì)象的有根系統(tǒng)發(fā)育對(duì)象Phylogeny[Rooted]中。為了生成系統(tǒng)發(fā)育樹(shù)??，我們將使用q2-phylogeny插件中的align-to-tree-mafft-fasttree工作流程。
首先，工作流程使用mafft程序執(zhí)行對(duì)FeatureData[Sequence]中的序列進(jìn)行多序列比對(duì)，以創(chuàng)建QIIME 2對(duì)象FeatureData[AlignedSequence]。接下來(lái)，流程屏蔽（mask或過(guò)濾）對(duì)齊的的高度可變區(qū)(高變區(qū))，這些位置通常被認(rèn)為會(huì)增加系統(tǒng)發(fā)育樹(shù)的噪聲。隨后，流程應(yīng)用FastTree基于過(guò)濾后的比對(duì)結(jié)果生成系統(tǒng)發(fā)育樹(shù)。FastTree程序創(chuàng)建的是一個(gè)無(wú)根樹(shù)，因此在本節(jié)的最后一步中，應(yīng)用根中點(diǎn)法將樹(shù)的根放置在無(wú)根樹(shù)中最長(zhǎng)端到端距離的中點(diǎn)，從而形成有根樹(shù)。
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](特征表[頻率])抽平到用戶(hù)指定的測(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距離（群落差異的定性度量，即只考慮種類(lèi)，不考慮豐度）
- 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)行無(wú)放回抽樣，從而使得結(jié)果表中的每個(gè)樣本的總計(jì)數(shù)為500。如果任何樣本的總計(jì)數(shù)小于該值，那么這些樣本將從多樣性分析中刪除。選擇這個(gè)值很棘手。建議通過(guò)查看上面創(chuàng)建的表table.qzv文件中呈現(xiàn)的信息并選擇一個(gè)盡可能高的值（因此每個(gè)樣本保留更多的序列）同時(shí)盡可能少地排除樣本來(lái)進(jìn)行選擇。
下面多樣性分析，需要基于重采樣/抽平(rarefaction)標(biāo)準(zhǔn)化的特征表，標(biāo)準(zhǔn)化采用無(wú)放回重抽樣至序列一致，如何設(shè)計(jì)樣品重采樣深度參數(shù)
--p-sampling-depth呢？
如是數(shù)據(jù)量都很大，選最小的即可。如果有個(gè)別數(shù)據(jù)量非常小，去除最小值再選最小值。比如此分析最小值為917，我們選擇1109深度重采樣，即保留了大部分樣品用于分析，又去除了數(shù)據(jù)量過(guò)低的異常值。本示例為近10年前測(cè)序技術(shù)的通量水平，454測(cè)序時(shí)代抽平至1000條即可，現(xiàn)在看來(lái)數(shù)據(jù)量很小。目錄一般采用HiSeq2500或NovaSeq6000的 PE250模式測(cè)序，數(shù)據(jù)量都非常大，通?？梢圆捎?萬(wàn)或5萬(wàn)的標(biāo)準(zhǔn)抽平，仍可保留90%以上樣本。過(guò)低或過(guò)高一般結(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
######################################################
