Tax4Fun是一個R包，可根據16S rRNA資料集預測微生物群落的功能（Aßhauer et, 2015）。Tax4Fun基於最小16S rRNA序列相似性的最近鄰識別實現16S rRNA基因序列與原核生物基因組功能註釋的連線。OTU丰度表必須來自SILVA特定版本的註釋，隨後Tax4Fun將OTU物種丰度轉換為功能配置檔案，該過程分為三個步驟。
首先，將基於SILVA的16S rRNA圖譜轉換為原核生物KEGG的分類學圖譜，線性變換是透過預先計算的關聯矩陣實現的。然後，根據NCBI基因組註釋獲得的16S rRNA複製數對KEGG生物的估計丰度進行歸一化處理。最後，利用歸一化的分類學丰度，將預測得到的KEGG生物的功能圖譜線性組合，預測微生物群落的功能圖譜。
本篇使用一個簡單示例，介紹使用R的Tax4Fun包，基於16S物種丰度資料預測細菌群落中代謝通路或功能基因的丰度，並在得到初步的預測結果之後，繼續使用R或STAMP分別對預測結果進行簡單的統計檢驗及作圖展示，以探索群落功能的方法。
瞭解Tax4Fun更多詳情請參見官網（Tax4Fun包也需在其官網中下載獲得）：http://tax4fun.gobics.de/
下文所有示例資料、執行結果檔案及R程式碼等，可在百度盤獲取（提取碼：si26）：
https://pan.baidu.com/s/1GVIBIsxIV5G_q4C3zApSjA
示例檔案
示例資料集共有80個16S測序樣本，均來自土壤。因試驗需求，在土壤中添加了某化學物質，目的為探究該化學物質對土壤微生物群落的影響。這80個樣本中，40個為不新增化學物質的對照組（control組），40個為新增化學物質的處理組（treat組）。
此處期望基於16S資料，根據群落物種組成使用Tax4Fun預測群落的功能，並想得知經過試驗處理後的細菌群落在功能上與對照組相比發生了怎樣的改變。
檔案“data/otu_table.txt”為OTU丰度表格，該OTU丰度表中無註釋列，即所有OTU的物種分類是未知的，在後續將使用這些OTU的代表序列（見下文）鑑定這些OTU的物種分類。
檔案“data/group.txt”為樣本分組資訊，第一列（names）為各樣本名稱；第二列（group）為各樣本的分組資訊，即這些樣本分別屬於未新增化學物質的對照組（control組）或添加了化學物質的處理組（treat組）。
檔案“data/otu.fasta”中包含了OTU丰度表中各OTU的代表序列。在後續將使用OTU代表序列與SILVA資料庫中的參考物種序列進行比對，以鑑定這些OTU所屬的“界門綱目科屬種”分類單元。
獲得Tax4Fun輸入檔案
因Tax4Fun要求輸入的OTU丰度表中，各OTU必須來自SILVA資料庫的註釋；並且Tax4Fun目前所具有3個參考版本，分別對應了3個SILVA物種註釋庫版本，最好二者相互對應。因此，為了預測結果的準確性，首先需要保證我們的OTU註釋資訊來自於SILVA資料庫的註釋。若最初的OTU註釋資訊非SILVA註釋（如來自Greengenes、RDP等），此時需要使用SILVA資料庫重新註釋。SILVA註釋時，最好使用與所將使用Tax4Fun參考版本所對應的SILVA物種註釋資料集版本。
例如，以下分析示例中挑選了Tax4Fun的SILVA123版本進行細菌群落功能預測，因此此處最好使用SILVA資料庫中的SILVA123物種註釋庫對我們的OTU進行物種註釋。
使用qiime為OTU添加註釋（SILVA123物種註釋）
對於非SILVA註釋的OTU註釋資訊，可以使用這些OTU的代表序列重新進行註釋。這裡使用到無註釋資訊的OTU丰度表（otu_table.txt）以及包含各OTU的代表序列檔案（otu.fasta），使用16S資料常用分析工具qiime來完成。
使用各OTU的代表序列，將它們與SILVA123資料庫中的代表物種序列進行比對，找到每個OTU序列的最相似序列，以期得到每個OTU所對應的物種分類單元（OTU註釋資訊）。同時將OTU丰度錶轉化為便於qiime識別的biom樣式，並使用qiime將各OTU註釋資訊新增至OTU表的最後一列。最終得到具有SILVA123物種註釋資訊的OTU丰度表格。
##注：以下步驟均在 shell 命令列中完成，執行程式或指令碼均來自 qiime
#在 otu_table.tsv 開頭新增一行“# Constructed from biom file”，以便將 otu_table.tsv 轉為 qiime 可識別樣式
cp otu_table.txt otu_table.tsv
sed -i '1i\# Constructed from biom file' otu_table.tsv
#otu_table.tsv 轉換為 otu_table.biom
biom convert -i otu_table.tsv -o otu_table.biom –table-type="OTU table" –to-json
#OTU 註釋，輸出結果 otu_table_tax_assignments.txt 即為註釋檔案
assign_taxonomy.py -i otu.fasta -r SILVA_123_SSURef_Nr99_tax_silva.fasta -t SILVA_123_SSURef_Nr99_tax_silva.tax -o ./
#新增 otu 註釋資訊至 biom 格式的 otu 表（otu_table.biom ）的最後一列，並將列名稱命名為 taxonomy
biom add-metadata -i otu_table.biom –observation-metadata-fp otu_table_tax_assignments.txt -o otu_table.silva.biom –sc-separated taxonomy –observation-header OTUID,taxonomy
#otu_table.silva.biom 轉換為 otu_table.silva.tsv
biom convert -i otu_table.silva.biom -o otu_table.silva.tsv –to-tsv –header-key taxonomy
網盤附件“qiime_result/otu_table.silva.tsv”為最後獲得的OTU丰度表，其第一行為“# Constructed from biom file”；最後一列為OTU註釋資訊，且註釋來源為SILVA資料庫；其餘行列即為各OTU名稱、各樣本名稱、以及各OTU在各樣本中的丰度分佈。
Tax4Fun預測細菌群落功能（預設流程）
以下為Tax4Fun進行功能預測的預設流程，功能預測結果將得到KEGG功能基因及代謝通路的丰度，並進行了KEGG功能基因及代謝註釋。
輸入了OTU丰度表，包含了每個測序樣本中物種型別組成及其丰度資訊。Tax4Fun會根據這些資訊，在其資料庫中找到對應的物種類別，根據16S複製數將物種資料標準化後，統計所對應物種類別其本身所具有的功能基因即所含代謝通路的數量，預測細菌群落中所具有的功能基因丰度以及代謝通路丰度。因此可知，OTU註釋結果中，若註釋資訊越詳細（分類水平越細，越接近種水平），則理論上功能註釋結果的準確度越高。
library(Tax4Fun)
##Tax4Fun 預設流程
#此處使用“otu_table.silva.tsv”作為輸入檔案
QIIMESingleData <- importQIIMEData('otu_table.silva.tsv')
folderReferenceData <- 'SILVA123'
#KEGG 功能基因（KO 第四級）丰度預測
Tax4FunOutput <- Tax4Fun(QIIMESingleData, folderReferenceData, fctProfiling = TRUE, refProfile = 'UProC', shortReadMode = TRUE, normCopyNo = TRUE)
tax4fun_gene <- as.data.frame(t(Tax4FunOutput$Tax4FunProfile))
gene <- rownames(tax4fun_gene)
write.table(cbind(gene, tax4fun_gene), 'tax4fun.gene.txt', row.names = FALSE, sep = '\t', quote = FALSE)
#KEGG 代謝通路（KO 第三級）丰度預測
Tax4FunOutput <- Tax4Fun(QIIMESingleData, folderReferenceData, fctProfiling = FALSE, refProfile = 'UProC', shortReadMode = TRUE, normCopyNo = TRUE)
tax4fun_pathway <- as.data.frame(t(Tax4FunOutput$Tax4FunProfile))
pathway <- rownames(tax4fun_pathway)
write.table(cbind(pathway, tax4fun_pathway), 'tax4fun.pathway.txt', row.names = FALSE, sep = '\t', quote = FALSE)
兩個流程裡面，一個重要引數“fctProfiling”可分別設定為“TRUE”和“FALSE”。當“fctProfiling = TRUE”時，將輸出功能基因丰度預測結果；當“fctProfiling = FALSE”時，將輸出代謝通路預測結果。其他引數一般情況下無需更改，預設即可。
對於功能基因丰度預測結果“tax4fun.gene.txt”（下圖1），第1列為各功能基因KO ID（可根據此ID在KEGG資料庫中檢索以檢視其詳細功能資訊等）、名稱及其對應的酶型別（酶命名規則可參考https://enzyme.expasy.org/enzyme_details.html），第2列及之後為各樣本中各功能基因的相對丰度。
對於代謝通路預測結果“tax4fun.pathway.txt”（下圖2），第1列為各KEGG代謝通路KO ID（可根據此ID在KEGG資料庫中檢索以檢視其詳細功能資訊等）及其名稱，第2列及之後為各樣本中各代謝通路的相對丰度。
以上即透過Tax4Fun功能預測得到了一種類似宏基因組測序分析中的KEGG代謝註釋結果。若其中有我們感興趣的代謝通路或者功能基因，或者想著重具體地瞭解某些代謝通路或者功能基因，可以進入KEGG官方網站（https://www.genome.jp/kegg/pathway.html），在網站中輸入這些通路或者基因的ID進行查詢，獲知其詳情。
例如，查詢代謝通路“Glycolysis / Gluconeogenesis”瞭解詳情，即“00010”，過程和結果如下所示。
根據預測結果進行統計分析（一個簡單的示例）
以及透過一些統計分析，觀測所預測出的代謝通路或功能基因丰度在各分組中是否具有顯著差異，以便後續有針對性地查閱相關文獻等。
例如以下展示一個小示例。
為 KEGG代謝通路結果新增分類註釋
例如，將KEGG代謝通路對映到更高階的KO（KEGG Orthology）分類單元（如第一、二級），這樣，相較於一開始的輸出結果（只包含KEGG代謝通路名稱），添加了這些代謝通路所屬的更高階的分類資訊。
#為 KEGG 代謝通路對映 KO 分類
#“pathway.anno.txt”為已經整理好的 KEGG 註釋檔案（可見網盤附件）
kegg_anno <- read.delim('pathway.anno.txt', sep = '\t', colClasses = 'character', check.names = FALSE)
tax4fun_pathway$KO3_id <- t(data.frame(strsplit(rownames(tax4fun_pathway), ';')))[ ,1]
tax4fun_pathway$KO3_id <- gsub('ko', '', tax4fun_pathway$KO3_id)
tax4fun_pathway <- merge(tax4fun_pathway, kegg_anno, by = 'KO3_id')
write.table(tax4fun_pathway, 'tax4fun.pathway.anno.txt', row.names = FALSE, sep = '\t', quote = FALSE)
KO1_id：KO第一級分類的ID；
KO1：KO第一級分類的名稱；
KO2_id：KO第二級分類的ID；
KO2：KO第二級分類的名稱；
KO3_id：KO第三級分類的ID；
KO3：KO第三級分類的名稱，也即具體的KEGG代謝通路；
其餘列：預測得到的丰度資訊。
KO二級分類單元的組間差異檢驗
繼續關注這些不同的功能在兩個不同分組（處理組與對照組）的差異。
以下不妨以KO第二級分類為例展示，計算對映到的丰度，並根據分組資訊計算差異。
#在KO 第二級統計求和，並新增樣本分組資訊
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
names(group)[1] <- 'variable'
tax4fun_pathway <- tax4fun_pathway[c(group$variable, 'KO2_id')]
tax4fun_pathway <- reshape2::melt(tax4fun_pathway, id = 'KO2_id')
tax4fun_pathway <- doBy::summaryBy(value~variable+KO2_id, tax4fun_pathway, FUN = sum)
tax4fun_pathway <- merge(tax4fun_pathway, group, by = 'variable')
#計算均值（mean）±標準差（sd），或標準誤差（se）
se <- function(x) sd(x) / (length(x))^0.5
pathway_stat <- doBy::summaryBy(value.sum~group+KO2_id, tax4fun_pathway, FUN = c(mean, sd, se))
#添加註釋
kegg_anno_2 <- kegg_anno[!duplicated(kegg_anno$KO2), ][-c(5:7)]
pathway_stat <- merge(pathway_stat, kegg_anno_2, by = 'KO2_id', all.x = TRUE)
#write.table(pathway_stat, 'tax4fun.KO2.anno.txt', row.names = FALSE, sep = '\t', quote = FALSE)
#顯著性差異分析，例如這裡直接使用非引數的 wilcoxon 檢驗兩組間差異
KO2_id <- unique(tax4fun_pathway$KO2_id)
for (i in KO2_id) {
tax4fun_pathway_2_i <- subset(tax4fun_pathway, KO2_id == i)
test <- wilcox.test(value.sum~group, tax4fun_pathway_2_i)
line_t <- which(pathway_stat$KO2_id == i & pathway_stat$group == 'treat')
pathway_stat[line_t,'p_value'] <- test$p.value
if (test$p.value < 0.05 & test$p.value >= 0.01) {
pathway_stat[line_t,'sign'] <- '*'
}
if (test$p.value < 0.01 & test$p.value >= 0.001) {
pathway_stat[line_t,'sign'] <- '**'
}
if (test$p.value < 0.001) {
pathway_stat[line_t,'sign'] <- '***'
}
}
#write.table(pathway_stat, 'tax4fun.KO2.anno_stat.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')
“tax4fun.KO2.anno_stat.txt”結果中包含了各KO功能分類單元的名稱，在各組中丰度的均值、標準差、標準誤差，以及wilcoxon檢驗的顯著性p值等資訊。
ggplot2作圖展示功能預測及組間差異分析結果
接下來，使用ggplot2進行結果的視覺化，展示本次Tax4Fun功能預測及顯著性差異分析結果。
#使用 ggplot2 作圖
library(ggplot2)
#給 KO 名稱按丰度數值大小排個序
pathway_stat <- pathway_stat[order(pathway_stat$value.sum.mean), ]
pathway_stat$KO2 <- factor(pathway_stat$KO2, levels = unique(pathway_stat$KO2))
pathway_stat$KO1 <- factor(pathway_stat$KO1, levels = rev(unique(pathway_stat$KO1)))
#作圖
pathway_stat$value.sum.mean <- 100 * pathway_stat$value.sum.mean
pathway_stat$value.sum.sd <- 100 * pathway_stat$value.sum.sd
ko2_plot <- ggplot(pathway_stat, aes(x = KO2, y = value.sum.mean, fill = group)) +
geom_col(position = 'dodge', width = 0.8, colour = 'black', size = 0.05) + #“dodge 柱狀圖”樣式
scale_fill_manual(values = c('red', 'blue')) + #填充顏色
geom_errorbar(aes(ymin = value.sum.mean – value.sum.sd, ymax = value.sum.mean + value.sum.sd),
size = 0.05, width = 0.35, position = position_dodge(width = 0.8)) + #新增誤差線（均值±標準差）
geom_text(aes(label = sign, y = value.sum.mean + value.sum.sd + 0.5), size = 4, position = position_dodge(0.8)) + #新增顯著性標記“*”
facet_grid(KO1~., space = 'free', scale = 'free_y') + #分面圖，按 KO1 分類
coord_flip() + #橫、縱座標軸反轉
theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black'),
legend.title = element_blank(), legend.position = 'top') + #背景、圖例主題等設定
labs(x = 'KEGG Orthology (KO)', y = 'Relative Abundance (%)') #座標軸標題
ggsave('ko2_stat.pdf', ko2_plot, width = 8, height = 10)
ggsave('ko2_stat.png', ko2_plot, width = 8, height = 10)
參考文獻
Aßhauer, Kathrin P, Wemheuer B , Daniel R , et al. Tax4Fun: predicting functional profiles from metagenomic 16S rRNA data. Bioinformatics, 2015, 31(17):2882-2884.
