PICRUSt（Phylogenetic Investigation of Communities by Reconstruction of Unobserved States）是一種使用標記基因資料和參考基因組資料庫來預測宏基因組功能組成的計算方法，可基於16S rRNA圖譜預測微生物群落的基因組功能組成（Langille et al, 2013）。
PICRUSt使用一種擴充套件的祖先狀態重建演算法來預測存在哪些基因家族，然後結合基因家族來估計複合宏基因組，包含一個兩步過程。在最初的“基因含量推斷”步驟中，透過參考系統發育樹預先計算每個生物體的基因含量，這將重建一個基於16S的系統發育中每種生物的預測基因家族丰度表。隨後的“宏基因組推斷”步驟將對所有微生物分類群的所得基因含量預測與16S rRNA基因的複製數進行的校正，最終獲得整個微生物群落預期的基因家族丰度。
接下來簡介PICRUSt工具的使用。
考慮到線上版PICRUSt的網上教程貌似挺多的，但本地版PICRUSt的使用很少有教程提到；並且PICRUSt的輸入檔案（GreenGene註釋的OTU表）也需本地操作，為了連貫性，所以本篇只介紹本地版PICRUSt工具。
PICRUSt簡介及官方教程：http://picrust.github.io/picrust/index.html
下文所有示例資料、執行結果檔案及命令列等，可在百度盤獲取（提取碼：sech）：
https://pan.baidu.com/s/1VsCThQWPddCvng02Pdmv3A
PICRUSt安裝
conda安裝PICRUSt環境
PICRUSt已經打包在bioconda裡面了，可直接透過conda安裝環境。
#使用 bioconda 安裝 PICRUSt 環境
conda create -n picrust1 -c bioconda -c conda-forge picrust
#啟用使用
source activate picrust1
#用完記得退出
source deactivate picrust1
手動下載需要的配置庫檔案
PICRUSt根據16S物種資訊，預測群落功能。在這一過程中，需要讀取相關的資料庫內容資訊，執行預測。PICRUSt提供了相關的庫檔案，根據需要下載。
連結：http://picrust.github.io/picrust/picrust_precalculated_files.html#id1
下載的檔案是“.gz”壓縮格式的，但無需解壓。將下載好的檔案放入PICRUSt環境安裝位置的data路徑下。
#比方說這裡使用 conda 安裝的 picrust 環境
#那麼，data 目錄就位於 conda 中的 picrust 環境目錄下
#具體而言，比方說我的 conda 安裝在
#~/software/Miniconda3
#那麼本示例中，透過 conda 安裝的 picrust 環境的 data 目錄，即位於
#~/software/Miniconda3/envs/picrust1/lib/python2.7/site-packages/picrust/data/
#下載好的一系列的庫檔案，如“16S_13_5_precalculated.tab.gz”等，就移至這裡即可
PICRUSt預測16S群落資料
軟體安裝並下載好配置檔案後，就可以根據的16S物種資料預測群落功能了。
1、準備GreenGene註釋的OTU丰度表
對於輸入的OTU丰度表，必須為GreenGene資料庫註釋的。並且，PICRUSt識別的並不是常規OTU表中最後一列的註釋資訊，而是OTU的GreenGene ID。具體而言，就是這種形式。
如果不清楚，可以透過如下方法獲得。
*透過qiime程式實現OTU序列的GreenGene註釋
OTU表不是GreenGene註釋的，就需要使用GreenGene進行註釋。這裡就以16S群落分析常用程式qiime為例進行註釋。
示例檔案可見網盤附件“data”。示例資料集共有80個16S測序樣本，均來自土壤。因試驗需求，在土壤中添加了某化學物質，目的為探究該化學物質對土壤微生物群落的影響。這80個樣本中，40個為不新增化學物質的對照組（control組），40個為新增化學物質的處理組（treat組）。
檔案“data/otu_table.txt”為OTU丰度表格，該OTU丰度表中無註釋列，即所有OTU的物種分類是未知的，在後續將使用這些OTU的代表序列（見下文）鑑定這些OTU的物種分類。
檔案“data/otu.fasta”中包含了OTU丰度表中各OTU的代表序列。在後續我們將使用OTU代表序列與GreenGene資料庫中的參考物種序列進行比對，以鑑定這些OTU所屬的“界門綱目科屬種”分類單元。
首先根據獲得的OTU的代表序列，透過GreenGene註釋，並將註釋結果新增到不帶註釋但已聚類好的OTU丰度表中。隨後再根據註釋概要，對映GreenGene ID，並用GreenGene ID替換原OTU表的OTU ID。
##OTU 序列註釋
#（1）首先解壓 gg_13_5_otus.tar.gz，這個是 GreenGene 資料庫 13.5 版本
#（2）藉助 qiime 程式實現 OTU 序列註釋
#如果沒安裝 qiime，在 linux 中，透過 conda 安裝 qiime 環境
#conda create -n qiime1 qiime
#啟用使用
#source activate qiime1
#用完記得退出
#source deactivate qiime1
#啟用 qiime 環境
#註釋完成後，先不著急退出 qiime，後面的某些命令還需用它，最後再退
source activate qiime1
#qiime 的 OTU 序列註釋過程
assign_taxonomy.py -i otu.fasta -r gg_13_5_otus/rep_set/97_otus.fasta -t gg_13_5_otus/taxonomy/97_otu_taxonomy.txt –uclust_max_accepts 1 -o gg97
#根據具體的註釋細節，對應 greengene id（這裡用了一個手寫的 R 指令碼實現轉換，在網盤附件中）
grep 'H' gg97/otu_tax_assignments.log > gg97/otu_tax_assignments2.log
Rscript replace_id.r gg97/otu_tax_assignments2.log otu_table.txt otu_table2.txt
#轉為 biom 格式
sed -i '1i\# Constructed from biom file' otu_table2.txt
biom convert -i otu_table2.txt -o otu_table2.biom –table-type="OTU table" –to-json
qiime的輸出可參考網盤附件“qiime_output”，“otu_table2.txt”即為最後獲得的以GreenGene ID為註釋的純文字的OTU表；“qiime_output/otu_table2.biom”為轉換為biom格式的OTU表，用於程式識別。
2、可選校正16S複製數（推薦）
對於很多細菌而言，一個個體可能包含多條16S（多複製16S）。PICRUSt提供了校正16S複製數的方法。
#啟用 picrust 環境
#已經激活了 qiime，再啟用 picrust，功能不衝突
source activate picrust1
#校正 16S 複製數
normalize_by_copy_number.py -i otu_table2.biom -o normalized_otus.biom
#轉為 txt 表格
biom convert -i normalized_otus.biom -o normalized_otus.txt –to-tsv
可參考網盤附件“picrust_output/normalized_otus.txt”，即為校正16S複製數後的OTU丰度表；“picrust_output/normalized_otus.biom”為轉換為biom格式的OTU表，用於程式識別。
3、群落功能預測
接下來，就使用校正16S複製數（物種丰度）後的OTU資料，執行功能預測。
#預測 KEGG 功能
predict_metagenomes.py -i normalized_otus.biom -o ko.biom -a nsti_ko.txt
#轉為 txt 表格（分別為 KO 層級和 KO 功能描述）
biom convert -i ko.biom -o ko_pathways.txt –to-tsv –header-key KEGG_Pathways
biom convert -i ko.biom -o ko_description.txt –to-tsv –header-key KEGG_Description
#預測 COG 功能
predict_metagenomes.py -i normalized_otus.biom -o cog.biom -a nsti_cog.txt –type_of_prediction cog
#轉為 txt 表格（分別為 COG 層級和 COG 功能描述）
biom convert -i cog.biom -o cog_category.txt –to-tsv –header-key COG_Category
biom convert -i cog.biom -o cog_description.txt –to-tsv –header-key COG_Description
#預測 Rfam 功能
predict_metagenomes.py -i normalized_otus.biom -o rfam.biom -a nsti_rfam.txt –type_of_prediction rfam
#轉為 txt 表格
biom convert -i rfam.biom -o rfam_description.txt –to-tsv –header-key RFAM_Description
上述結果均可在網盤附件“picrust_output”中找到。
PICRUSt提供了對KEGG、COG、Rfam等資料庫功能的對映。預設情況下，僅預測KEGG；其它資料庫功能的預測，透過“–type_of_prediction”引數指定類別。
PICRUSt預設輸入輸出檔案均為biom格式，可再轉換為純文字樣式，檔案結構類似OTU丰度表。以COG預測結果為例，第一列為預測所得COG功能ID，最後一列記錄了該COG所示分類或功能描述，中間列為該功能的丰度資訊。
獲得群落功能丰度表後，就可以按照OTU丰度表的統計分析方法，去執行類似的分析了。這點可以找一些文獻作參考，看別人是怎樣做的。例如，首先計算特定功能丰度在組間的顯著性，獲得組間差異顯著的功能，然後再從資料庫官網上查詢該功能的細節，解釋生物學現象等。
4、PICRUSt的一些其它實用功能
除了獲得整個群落水平的功能丰度外，PICRUSt還能提供在不同功能分類水平上的統計，以及單個OTU對功能的貢獻。
#分解 KO，比方說在 KO 第 2 級或第 3 級功能上統計
categorize_by_function.py -f -i ko.biom -c KEGG_Pathways -l 3 -o ko_L3.txt
categorize_by_function.py -f -i ko.biom -c KEGG_Pathways -l 2 -o ko_L2.txt
#分解 COG，比方說在 COG 第 2 級功能上統計
categorize_by_function.py -f -i cog.biom -c COG_Category -l 2 -o cog_L2.txt
#OTU 對功能的貢獻，需指定特定 KO 功能分類 id 或 COG 功能分類 id
metagenome_contributions.py -i normalized_otus.biom -l K00001,K00002,K00004 -o ko_metagenome_contributions.txt
metagenome_contributions.py -i normalized_otus.biom -l COG0001,COG0002 -t cog -o cog_metagenome_contributions.txt
#用完記得退出 qiime 或 picrust 環境哦
source deactivate picrust1
source deactivate qiime
上述結果均可在網盤附件“picrust_output”中找到。
biom檔案轉換為純文字檔案才可方便檢視，以下繼續以COG的預測結果為例展示。總之PICRUSt使用起來挺方便的，特別是它自己就含有直接在各功能分類層級上的統計功能，這樣倒是可以省去再使用R等工具去作分類合併了。如果期望關注某些OTU是否對群落功能是重要的，貢獻度表可以提供參考。
隨後可能的統計分析方法，參考實際的文獻中的方法來就可以了。
不妨作個圖觀察下組間功能分類的差異吧，還是以COG為例。
#R 語言操作
#讀取 COG 二級分類統計
#把“cog_L2.txt”第一行的“# Constructed from biom file”刪除再讀取
cog <- read.delim('cog_L2.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
cog <- cog[-ncol(cog)]
#整理分類註釋
cog$category <- apply(cog[1], 1, function(x) substr(x, 2, 2))
names(cog)[1] <- 'description'
cog <- reshape2::melt(cog, id = c('category', 'description'))
#合併樣本分組檔案
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
names(group)[1] <- 'variable'
cog <- merge(cog, group, by = 'variable')
#ggplot2 作圖，細節就不調整了
library(ggplot2)
ggplot(cog, aes(x = category, y = value)) +
geom_boxplot(aes(color = group))
參考文獻
Langille M G I, Zaneveld J, Caporaso J G, et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nature Biotechnology, 2013, 31(9):814-821.
