中文整合包_MIMOSA2： 基于微生物组和代谢组数据的整合分析
MIMOSA2：基于微生物組和代謝組數據的整合分析
MIMOSA2 升級自MIMOSA1。是 Borenstein 實驗室(http://borensteinlab.com/ , 專注宏基因組系統 生物學)最新開發的工具。用于微生物群落和代謝組的整合分析，尋找微生物和代謝產物之間的關系。
先前Borenstein bab開發過Fishtaco工具(http://borensteinlab.com/software_fishtaco.html)，用于微生物群落和功能聯合分析，[Fishtaco工具無論是擴增子測序數據還是宏基因組數據都可以使用，只是功能和微生物群落數據的對應關系不同，詳細參見我自己寫的Fishtaco教程，這也是目前唯一一份中文教程:參見Fishtaco](https://mp.weixin.qq.com/s/DBSFnTYGofRgEgGn7QGWjQ)。本次MIMOSA2用于整合微生物群落并在功能層面上耦合代謝數據，希望可以找到關鍵微生物和特定代謝產物的之間的關系。作者為此還開發了網頁版本的工具：https://elbo-spice.gs.washington.edu/shiny/MIMOSA2shiny/ ，可以先嘗鮮。
MIMOSA2有一個數據庫，包含物種及其對應的代謝能力信息，我們的擴增子測序數據可以通過不同的方式映射至數據庫。讓我欣喜的是無論是基于去噪方式的ASV還是傳統聚類方式的OTU，都可以很好地進行分析。這也是隨著新的聚類方式發展的新式工具之一。
MIMOSA2通過微生物群落特征和代謝組特征關聯，回答以下幾個問題：
代謝組數據是否和微生物組數據緊密相關；
微生物群落差異是否可以解釋代謝差異；
何種微生物或者基因造成了代謝差異，或者貢獻了代謝差異。
MIMOSA2的工作原理
工作原理：首先得到微生物組數據，使用參考數據庫構建代謝模型；計算不同微生物對不同代謝物的潛在代謝潛能指標(CMP)；對每種代謝產物全部潛在貢獻微生物CMP值進行匯總并做回歸分析，評估微生物群落對指定代謝物的預測能力；最后篩選對特定代謝物影響加大的微生物作為關鍵微生物。
MIMOSA2分析的主要步驟
構建群落和代謝模型，這里包括群落中物種對代謝通路的貢獻程度；
使用代謝模型計算每種微生物及其代謝得，評估在每個樣本匯總每種微生物和代謝產物的影響。
計算群落水平的代謝潛能得分，使用回歸模型評估在不同樣本中是否差異。
可視化特征微生物對特定代謝的影響，并尋找關鍵微生物。
作者做了一份完整的介紹，參見網頁：https://borenstein-lab.github.io/MIMOSA2shiny/analysis_description.html
上圖展示了MIMOSA2分析所有可能的組合分析形式：使用微生物組和宏基因組數據作為輸入，通過MIMOSA2可選構建兩個模型，分別為：AGORA和EMBL；本次我為大家演示擴增子數據和代謝組數據如何構建模型并進行后續分析。從圖中我們可以看到基于Greengenes數據庫和SILVA數據庫的物種注釋結果，還有目前最為流行的非聚類方式得到的ASV(amplicon sequence variant)都可以用于分析。這里我們可以看到需要兩個數據庫：AGORAh和EMBL數據庫。所以首先我們來下載數據庫。
下面我們下載數據庫，基于非聚類方式ASV的兩種模型數據庫和基于傳統聚類(OTU)的兩種模型數據庫，注意這里基于OTU的兩個模型適用于greengene和Silva數據庫注釋結果。download_reference_data為數據庫下載命令，網站 https://borenstein-lab.github.io/MIMOSA2shiny/downloads.html# 可以手動下載數據庫，但是好像我不能下載，顯示拒絕訪問。大家還是用下面命令好些。數據庫默認下載到當前文件夾中，所以如果后使用多的話指定一個固定文件夾。四個數據庫下載下來大小一共為500多M。
軟件部署
安裝R包
```{R warning=TRUE, include=FALSE}
rm(list=ls())
#安裝包 ? 這里我安裝完成之后就將其注釋掉
library(devtools)
install_github(“borenstein-lab/mimosa2”)
#如果安裝失敗，下載github包，進行本本地安裝：
install.packages(“C:/Users/liulanlan/Desktop/mimosa2-master”, repos = NULL, type = “source”)
載入依賴關系```{R}
suppressWarnings(suppressMessages(library(ggcorrplot)))
suppressWarnings(suppressMessages(library(igraph)))
suppressWarnings(suppressMessages(library(psych)))
suppressWarnings(suppressMessages(library(network)))
suppressWarnings(suppressMessages(library(ggplot2)))
suppressWarnings(suppressMessages(library(sna)))
suppressWarnings(suppressMessages(library(ergm)))
suppressWarnings(suppressMessages(library(ggrepel)))
suppressWarnings(suppressMessages(library("mimosa")))
suppressWarnings(suppressMessages(library(data.table)))
下載參考數據庫
```{R include=FALSE}
基于非聚類方法
download_reference_data(seq_db = “Sequence variants (ASVs)”, target_db = “AGORA genomes and models”)
download_reference_data(seq_db = “Sequence variants (ASVs)”, target_db = “RefSeq/EMBL_GE Ms genomes and models”)
#基于參考數據庫方法
download_reference_data(seq_db = “Greengenes 13_5 or 13_8 OTUs”, target_db = “RefSeq/EMBL_GEMs genomes and models”)
download_reference_data(seq_db = “Greengenes 13_5 or 13_8 OTUs”, target_db = “AGORA genomes and models”)
### 我們來看看MIMOSA2究竟做了什么？
A. M為模擬的物質，M濃度的變化和物種1-3之間的關系如下圖所示。
B. A-F為不同樣本，物種1和2對于M物質具有相同的變化，這表明了這兩種物質促進M物質形成，物種3消耗M物質。比較群落水平的代謝潛力(Compare total community-level metabolic potential, CMP)是最重要的指標，代表每個物種對該物質的影響得分。
C. 表示CMP進行回歸并評估對給物質影響的顯著性。這里給出R為評估的全部物種對該物質的解釋度。
D. 給出了三個物種影響程度，這里最大的`Taxon 3`也就是所謂的支配物種，對M物質具有支配作用。
### 運行計算
這里一共有四種模式的運算：
#### 模式1：基于Greengenes OTU的結果和EMBL模型
第一種模式：使用基于聚類的結果，也就是OTU進行分析，我們在這里選擇RefSeq/EMBL_GEMs genomes and models模型，其對應的數據庫再之前下載得到為OTU_EMBL，這里我們指定數據庫的相對路徑，并指定到其中的data中。
參數文件：config_example1.txt，方便理解查查看，這也是我們唯一的輸入，因為全部的參數都在這里指定好了。這里值得注意的參數是vsearch_path，我們需要下載vsearch并安裝，這里的安裝無論是再linux還是win下直接賦予可執行權限并放入環境變量中即可。本例子我運行在win環境，這里必須寫全稱：vsearch.exe，在linux下只需要寫vsearch就好了。
參數文件(復制保存為txt即可調用)：
```{R eval=FALSE, include=FALSE}
file1 test_gg.txt #gg參考數據庫聚類結果
file2 test_mets.txt#代謝物質結果，使用KEGG編號
ref_choices RefSeq/EMBL_GEMs genomes and models#使用模型
file1_type Greengenes 13_5 or 13_8 OTUs#定義輸入微生物數據類型
simThreshold 0.99#相似度
data_prefix ./database/OTU_EMBL/data/#定義參考數據庫位置
vsearch_path vsearch.exe#定義vsearch位置，此時在win下進行，所以指定為exe
rank_based T
metType KEGG Compound IDs#指定代謝產物ID類型
signifThreshold 1
運行第一種模式
# 下面我們選擇一個例子進行運行，全部輸入文件和參考數據庫都備注再這里：config_example1.txt#---運行-首先測試第一種模式：就gg數據庫聚類結果使用RefSeq/EMBL_GEMs genomes and models模型分析
mimosa_results = run_mimosa2("./config_example1.txt", make_plots = T, save_plots = T)
# unclass(mimosa_results)
模式2：基于Greengenes OTU和AGORA模型
第二種模式:使用基于聚類的結果，也就是OTU進行分析，我們在這里選擇AGORA genomes and models模型，其對應的數據庫再之前下載得到為OTU_AGORA，這里我們指定干數據庫的相對路徑，并指定到其中的data中，
參數文件：./test4.txt
```{R eval=FALSE, include=FALSE}
file1 ? ?test_gg.txt
file2 ? ?test_mets.txt
ref_choices ? ?RefSeq/EMBL_GEMs genomes and models
file1_type ? ?Greengenes 13_5 or 13_8 OTUs
simThreshold ? ?0.99
data_prefix ? ?./database/OTU_AGORA/data/
vsearch_path ? ?vsearch.exe
rank_based ? ?T
metType ? ?KEGG Compound IDs
signifThreshold ? ?1
#---運行=第二個測試基于OTU結果的RefSeq/EMBL_GEMs genomes and models模型分析
mimosa_results = run_mimosa2("./test4.txt", make_plots = T, save_plots = T)
模式3：基于ASV和EMBL模型
第三種模式: 使用基于非聚類的結果，也就是ASV進行分析，file1的文件類型需要更改為Sequence variants (ASVs)，我們在這里選擇EMBL模型，其對應的數據庫再之前下載得到為ASV_EMBL，這里我們指定干數據庫的相對路徑，并指定到其中的data中，
參數文件：test3.txt
```{R eval=FALSE, include=FALSE}
file1 ? ?test_seqs.txt
file2 ? ?test_mets.txt
ref_choices ? ?EMBL_GEMs genomes and models
file1_type ? ?Sequence variants (ASVs)
simThreshold ? ?0.99
data_prefix ? ?./ASV_EMBL/data/
?vsearch_path ? ?vsearch.exe
rank_based ? ?T
metType ? ?KEGG Compound IDs
signifThreshold ? ?1
#---運行=第二個測試基于非聚類結果
mimosa_results = run_mimosa2("./test3.txt", make_plots = T, save_plots = T)
模式4：基于ASV和AGORA模型
第四種模式：使用基于非聚類的結果，也就是ASV進行分析，file1的文件類型需要更改為Sequence variants (ASVs)，我們在這里選擇AGORA genomes and models模型，其對應的數據庫再之前下載得到為ASV_AGORA，這里我們指定干數據庫的相對路徑，并指定到其中的data中，
參數文件：test4.txt
```{R eval=FALSE, include=FALSE}
file1 ? ?test_seqs.txt
file2 ? ?test_mets.txt
ref_choices ? ?AGORA genomes and models
file1_type ? ?Sequence variants (ASVs)
simThreshold ? ?0.99
data_prefix ? ?./database/ASV_AGORA/data/
vsearch_path ? ?vsearch.exe
rank_based ? ?T
metType ? ?KEGG Compound IDs
signifThreshold ? ?1
```{R}
#---運行=第二個測試基于非聚類結果的AGORA genomes and models模型分析
library(RColorBrewer)#調色板調用包
library("cowplot")
library(data.table)
# 主要是build_metabolic_model函數指定文件夾錯誤，所以我修改了該函數兩處文件夾的位置
source("./build_metabolic_model-wt.R")
#載入修改后的主函數run_mimosa2-wt
source("./run_mimosa2-wt.R")
#當我調整之后，進行分析已經沒有問題了
mimosa_results = run_mimosa2("./test2.txt", make_plots = T, save_plots = T)
結果說明
表格1: 微生物對代謝物的貢獻表格
最重要的一列：varshare，代表了該模型匯總為微生物可以解釋代謝物多少差異。僅僅包括貢獻p值小于0.1代謝物。也就是說這些代謝物受到那些微生物影響，或者那些微生物影響了這些物質變化。找出對物質變化影響比較大的微生物，作為后續研究方向。
表格2：原始物種表格匹配到數據庫后的新的物種表格
我們無論是基于ASV表格還是參考數據庫聚類的OTU表，那種模型，都會重新匹配模擬到新的物種，這個文件展示相關信息(我們可以看到雖然參加分析的OTU表樣品數量很多，但是僅展示了和代謝組共有的樣本)
表格3：CMP值矩陣
CMP值矩陣，展示了每個樣本中每種物種對每個物質的CMP評估值
表格4：模型統計摘要
模型統計總體簡要統計：參與分析的樣本數量，物種數量，參與模型的物種數量，代謝物質數量，參與模型的代謝物質數量，計算得到CMP值的代謝物質數量，匹配模型的代謝物質數量，顯著的代謝物質數量···
表格5：輸入參數文件，這里包含了我么的輸入文件，模型選擇，數據庫指定等信息
圖片文件展示
這里兩種物質，對應兩個結果。第一列圖片展示的是多個樣本涉及特定代謝物質的CMP值加和做的回歸，R值展示為加粗數字，數值越大代表微生物群落的預測能力越準確。第二列是對應的微生物的貢獻柱狀圖，柱狀圖右邊的是促進該物質生成，左邊的是消耗該物質。
手動輸出選擇：以上兩種方式我們得到的默認圖表如果是合并圖例，如果想單獨展示，這里結果中會包含圖片文件我們可以手動輸出，我在這里提供手動輸出的文件：
#手動輸出CMP擬合代謝物的圖表p = mimosa_results$CMPplots[[1]]
p+ theme_bw()#手動導出微生物對代謝物的貢獻圖
p = mimosa_results$metContribPlots[[1]]
p+ theme_bw()
輸出網絡：
這里我們可以看到基于微生物和代謝物之間建立了一個網絡，這是一個有向網絡，varshare值在這里類似我們傳統網絡中的R值。構建網絡的信息位于varshare值矩陣，varshare值作為邊的權重文件，我們可以添加對微生物群落對代謝物預測能力的R值。
首先我們來構造邊和節點文件，然后輸出，可以選擇的其他軟件進行網絡繪制：
# 構造邊文件networ = as.data.frame(mimosa_results$varShares)
# 構造一個向量用于承接正負關系
aa = rep(1,length(networ$VarShare))
for (i in 1:length(networ$VarShare)) {
if (networ$VarShare[i]> 0) {
aa[i] = 1
}else{
aa[i] = -1
}
}
# aa
# 組合成邊文件，添加了varshare值和方向。
edge = data.frame(compound = networ$compound,tax =networ$Species,value = networ$VarShare,direct= "directed",N_P = aa)
# head(edge)
#構造節點文件，代謝物和微生物共同組成，對分泌物的預測預測效果R值標記為一列
node = data.frame(ID= c(unique(networ$compound), unique(networ$Species) ),R = c(unique(networ$Rsq),rep("NA",length(unique(networ$Species)))))
# head(node)
write.csv(edge,"./Gephi_edge.csv")
write.csv(node,"./Gephi_node.csv")
R 中網絡可視化(可選)
#構造網絡文件G #轉化為鄰接矩陣，attr：填充權重
occor.r
#構造網絡文件，network包提供
g # (summary(g))
#做沒有權重的鄰接矩陣
m #設置可視化layout
plotcord #添加表頭
colnames(plotcord) = c("X1", "X2")
#添加節點名稱
plotcord$elements #制作邊文件
edglist edglist = as.data.frame(edglist)
# head(edglist)
#添加邊的權重
set.edge.value(g,"weigt",occor.r)
edglist$weight = as.numeric(network::get.edge.attribute(g,"weigt"))
# dim(edglist)
#添加其他信息
edges # head(edges)
edges$weight = as.numeric(network::get.edge.attribute(g,"weigt"))
##這里將邊權重根據正負分為兩類
aaa = rep("a",length(edges$weight))
for (i in 1:length(edges$weight)) {
if (edges$weight[i]> 0) {
aaa[i] = "+"
}
if (edges$weight[i]< 0) {
aaa[i] = "-"
}
}
#添加到edges中
edges$wei_label = as.factor(aaa)
colnames(edges) # head(edges)
##plotcord這個表格添加注釋
row.names(plotcord) = plotcord$elements
row.names(node) = node$ID
#合并之前制作的節點文件
index = merge(plotcord,node,by = "row.names",all = T)
# dim(index)
# head(index)
index$R = as.numeric(as.character(index$R))
#這里我為了突出代謝物，將全部微生物節點大小定義為代謝物的十分之一
index$R [is.na(as.numeric(as.character(index$R)))] = min(as.numeric(as.character(index$R)),na.rm = TRUE)/10
plotcord = index
#網絡出圖
pnet data = edges, size = 0.5) +
geom_point(aes(X1, X2,size=R),pch = 21,color = "black",fill = "red", data = plotcord) + scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
geom_text_repel(aes(X1, X2,label=elements),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
ggsave("./network.pdf", pnet, width = 12, height =8)
撰文：文濤 南京農業大學
責編：劉永鑫 中科院遺傳發育所
猜你喜歡
10000+：菌群分析?寶寶與貓狗?梅毒狂想曲 提DNA發Nature?Cell專刊?腸道指揮大腦
系列教程：微生物組入門 Biostar 微生物組 ?宏基因組
專業技能：學術圖表?高分文章?生信寶典 不可或缺的人
一文讀懂：宏基因組 寄生蟲益處 進化樹
必備技能：提問 搜索 ?Endnote
文獻閱讀 熱心腸 SemanticScholar Geenmedical
擴增子分析：圖表解讀 分析流程 統計繪圖
16S功能預測 ? PICRUSt ?FAPROTAX ?Bugbase Tax4Fun
在線工具：16S預測培養基 生信繪圖
科研經驗：云筆記 ?云協作 公眾號
編程模板:?Shell ?R Perl
生物科普:??腸道細菌?人體上的生命?生命大躍進 ?細胞暗戰 人體奧秘 ?
寫在后面
為鼓勵讀者交流、快速解決科研困難，我們建立了“宏基因組”專業討論群，目前己有國內外5000+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，并掃碼加主編好友帶你入群，務必備注“姓名-單位-研究方向-職稱/年級”。PI請明示身份，另有海內外微生物相關PI群供大佬合作交流。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍未解決群內討論，問題不私聊，幫助同行。
學習16S擴增子、宏基因組科研思路和分析實戰，關注“宏基因組”
總結
以上是生活随笔為你收集整理的中文整合包_MIMOSA2： 基于微生物组和代谢组数据的整合分析的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 防火墙IPSec 虚拟专用网络配置[虚拟
- 下一篇： CentOS7 编译安装LVS 互为主备