FAPROTAX是一款適用于對環(huán)境樣本(如海洋、湖泊等)的生物地球化學(xué)循環(huán)過程(特別是碳、氫豫喧、氮、磷幢泼、硫等元素循環(huán))進(jìn)行功能注釋預(yù)測的軟件紧显。我們知道微生物與土壤養(yǎng)分的循環(huán)利用密切相關(guān),許多養(yǎng)分被植物利用過程中離不開微生物的作用缕棵,比如鸟妙,氮降解途徑中的亞硫酸鹽呼吸。
FAPROTAX原理
首先作者先根據(jù)文獻(xiàn)資料手動(dòng)構(gòu)建了聯(lián)系物種分類與功能注釋的FAPROTAX數(shù)據(jù)庫挥吵;后又編寫了聯(lián)系OTU分類表與FAPROTAX數(shù)據(jù)庫的python腳本重父;最后,只要將基于16S的OTU分類表通過python腳本就可以輸出微生物群落功能注釋預(yù)測結(jié)果忽匈。
另外房午,PICRUSt輸入的物種豐度表目前只能用Greengene數(shù)據(jù)庫進(jìn)行物種注釋,在這一點(diǎn)上FAPROTAX更靈活丹允,因?yàn)樗R別的是菌的屬郭厌、種的名稱,只要你注釋到的菌己被報(bào)導(dǎo)有這方面的功能雕蔽,無論注釋數(shù)據(jù)庫是Greengene折柠、Silva和RDP,OTU是有參考還是de novo的結(jié)果都可以識別批狐。
有關(guān)原理部分和使用部分可參考宏基因組—如何用FAPROTAX預(yù)測微生物群落功能扇售。網(wǎng)址:https://mp.weixin.qq.com/s/J8EwJD_PTDhqRaD7kXlK1A
但是出現(xiàn)一個(gè)問題上面這篇文章給出的FAPROTAX官網(wǎng)網(wǎng)址(https://www.zoology.ubc.ca/louca/FAPROTAX/lib/php/index.php?section=Instructions)似乎失效了前塔,無法下載該軟件。
FAPROTAX官網(wǎng)的最新網(wǎng)址:https://pages.uoregon.edu/slouca/LoucaLab/archive/FAPROTAX/lib/php/index.php
可能需要翻墻才能訪問承冰,沒法翻墻的小伙伴也別急华弓,關(guān)注公眾號:Amoy數(shù)據(jù)分析,回復(fù)FAPROTAX即可下載最新版本的軟件困乒。
可以看到關(guān)于FAPROTAX的介紹寂屏,原理,使用代碼娜搂,如何下載選項(xiàng)迁霎,點(diǎn)擊下載選項(xiàng)可以看最新版本為1.2.1,點(diǎn)擊下載百宇。
軟件的安裝和使用
將下載的壓縮文件上傳至服務(wù)器或虛擬機(jī)考廉。這款軟件所需要的環(huán)境為Python2.7版本。
#解壓縮該文件
unzip FAPROTAX_1.2.1.zip
cd FAPROTAX_1.2.1 | ls
#使用conda軟件創(chuàng)建名為envpython27的python=2.7環(huán)境恳谎,-n為創(chuàng)建環(huán)境的名稱芝此。
conda create -n envpython27 python=2.7 -y
#激活python=2.7環(huán)境
conda activate envpython27
#測試軟件能否正確運(yùn)行憋肖,如果可以運(yùn)行會(huì)出現(xiàn)參數(shù)信息因痛,報(bào)錯(cuò)說明不可用,再根據(jù)報(bào)錯(cuò)的信息解決岸更。
python /home/ubuntu/FAPROTAX_1.2.1/collapse_table.py
#使用:-i otu分類信息表:-o 輸出文件名稱:-g 自帶的數(shù)據(jù)庫名稱鸵膏,在FAPROTAX_1.2.1/目錄下:-r 中間過程文件,可以看到那些菌屬?zèng)]有被數(shù)據(jù)庫對比上怎炊,其它具體使用請參考官網(wǎng)的教程谭企。
python /home/ubuntu/FAPROTAX_1.2.1/collapse_table.py -i qiime1_otu.txt -o qii1_otu_tax_faprtax.txt -g /home/ubuntu/FAPROTAX_1.2.1/FAPROTAX.txt -r report.txt -v -d 'taxonomy'
輸入的otu分類表文件格式:案列中的qiime1_otu.txt文件
輸出文件qii2_otu_tax_faprtax.txt
使用R語言進(jìn)行可視化
以下操作在R語言中進(jìn)行,一共需要三個(gè)文件:上一步的輸出文件qii2_otu_tax_faprtax.txt评肆;分組文件design.txt债查;需要展示的功能列表文件IND.list。后兩個(gè)文件自己創(chuàng)建即可瓜挽。
#加載包和主題盹廷,我們做圖的最終目的是要發(fā)表的,統(tǒng)一的設(shè)置一個(gè)主題可以方便我們作圖久橙,圖形更加美觀
library(ggplot2)
main_theme = theme(panel.background=element_blank(), panel.grid=element_blank(),
axis.line.x=element_line(size=.5, colour="black"), axis.line.y=element_line(size=.5, colour="black"),
axis.ticks=element_line(color="black"), axis.text=element_text(color="black", size=),
legend.position="right", legend.background=element_blank(), legend.key=element_blank(), legend.text=
element_text(size=),
text=element_text(family="sans", size=))
alpha = read.table("qii2_otu_tax_faprtax.txt", header=T, row.names=1, sep="\t", comment.char="")
alpha =alpha[,-1]#數(shù)據(jù)第一列是NA俄占,所以要去除
alpha = t(alpha)
alpha = alpha/100
alpha =as.data.frame(alpha)
design = read.table("design.txt", header=T, row.names=1, sep="\t")#SampleID為樣品編號與qii2_otu_tax_faprtax.txt第一行一致:subspecies為分組信息,對照組和處理組:solitype代表兩個(gè)不同地點(diǎn)淆衷。
design$group=design$groupID
sub_design = subset(design, group %in% c("MSXC","MSXT","XWZC","XWZT"))
sub_design$group = factor(sub_design$group, levels=c("MSXC","MSXT","XWZC","XWZT"))
idx = rownames(sub_design) %in% rownames(alpha)
sub_design=sub_design[idx,]
# detach("package:ggpubr", unload=TRUE)
# library(dplyr)
list = read.table("IND.list", header=F, sep="\t")#我們需要挑出來展示的功能缸榄,我關(guān)注的點(diǎn)是元素循環(huán)這塊,好像還有病原菌相關(guān)的祝拯,可根據(jù)自己的需要挑選
sub_alpha=alpha[rownames(sub_design), as.vector(list$V1)]
sub_alpha$sampleID=rownames(sub_alpha)
melted = melt(sub_alpha, id.vars = "sampleID")
melted_all = merge(melted, sub_design[,c("subspecies","solitype")], by.x="sampleID", by.y = "row.names", all.x = T )
x = as.data.frame(melted_all %>% group_by(variable) %>%
summarise(mean = mean(value)))
x = x[order(x$mean, decreasing = T),]#排序
x = head(x, 10)#選取排名前10的進(jìn)行展示
melted_all = melted_all[melted_all$variable %in% x$variable,]
melted_all$variable = factor(melted_all$variable, levels = as.vector(x$variable))
melted_all$soiltype = factor(melted_all$solitype, levels=c("M", "X"))
p = ggplot(melted_all, aes(x=variable, y=value, color=subspecies)) +
geom_boxplot(position = "dodge", alpha=1, outlier.size=0.3, size=0.5, width=0.7, fill="transparent") +
#geom_jitter( position=position_jitter(0.17), size=0.3, alpha=0.7)+
main_theme +
facet_grid(soiltype ~ .)+#以soiltype 圖型分頁
theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))
p
ggsave(paste("faprotax_IND_boxplot_top9..", ".pdf", sep=""), p, width = 8, height = 5)
design.txt文件格式甚带。
除了使用柱狀圖展示,還可使用熱圖等其它方式展示欲低。
致謝
公眾號宏基因組———如何用FAPROTAX預(yù)測微生物群落功能
R語言代碼來自文章———Zhang J , Liu Y X , Zhang N , et al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice[J]. Nature Biotechnology, 2019.中Fig3中的代碼部分辕宏。