MIMOSA2：基于微生物組和代謝組數(shù)據(jù)的整合分析

MIMOSA2 升級自MIMOSA1。是 Borenstein 實驗室(http://borensteinlab.com/ , 專注宏基因組系統(tǒng) 生物學)最新開發(fā)的工具。用于微生物群落和代謝組的整合分析，尋找微生物和代謝產(chǎn)物之間的關(guān)系。

先前Borenstein bab開發(fā)過Fishtaco工具(http://borensteinlab.com/software_fishtaco.html)，用于微生物群落和功能聯(lián)合分析，[Fishtaco工具無論是擴增子測序數(shù)據(jù)還是宏基因組數(shù)據(jù)都可以使用，只是功能和微生物群落數(shù)據(jù)的對應(yīng)關(guān)系不同，詳細參見我自己寫的Fishtaco教程，這也是目前唯一一份中文教程:參見Fishtaco](https://mp.weixin.qq.com/s/DBSFnTYGofRgEgGn7QGWjQ)。本次MIMOSA2用于整合微生物群落并在功能層面上耦合代謝數(shù)據(jù)，希望可以找到關(guān)鍵微生物和特定代謝產(chǎn)物的之間的關(guān)系。作者為此還開發(fā)了網(wǎng)頁版本的工具：https://elbo-spice.gs.washington.edu/shiny/MIMOSA2shiny/ ，可以先嘗鮮。

MIMOSA2有一個數(shù)據(jù)庫，包含物種及其對應(yīng)的代謝能力信息，我們的擴增子測序數(shù)據(jù)可以通過不同的方式映射至數(shù)據(jù)庫。讓我欣喜的是無論是基于去噪方式的ASV還是傳統(tǒng)聚類方式的OTU，都可以很好地進行分析。這也是隨著新的聚類方式發(fā)展的新式工具之一。

MIMOSA2通過微生物群落特征和代謝組特征關(guān)聯(lián)，回答以下幾個問題：

代謝組數(shù)據(jù)是否和微生物組數(shù)據(jù)緊密相關(guān)；

微生物群落差異是否可以解釋代謝差異；

何種微生物或者基因造成了代謝差異，或者貢獻了代謝差異。

MIMOSA2的工作原理

工作原理：首先得到微生物組數(shù)據(jù)，使用參考數(shù)據(jù)庫構(gòu)建代謝模型；計算不同微生物對不同代謝物的潛在代謝潛能指標(CMP)；對每種代謝產(chǎn)物全部潛在貢獻微生物CMP值進行匯總并做回歸分析，評估微生物群落對指定代謝物的預測能力；最后篩選對特定代謝物影響加大的微生物作為關(guān)鍵微生物。

MIMOSA2分析的主要步驟

構(gòu)建群落和代謝模型，這里包括群落中物種對代謝通路的貢獻程度；

使用代謝模型計算每種微生物及其代謝得，評估在每個樣本匯總每種微生物和代謝產(chǎn)物的影響。

計算群落水平的代謝潛能得分，使用回歸模型評估在不同樣本中是否差異。

可視化特征微生物對特定代謝的影響，并尋找關(guān)鍵微生物。

作者做了一份完整的介紹，參見網(wǎng)頁：https://borenstein-lab.github.io/MIMOSA2shiny/analysis_description.html

上圖展示了MIMOSA2分析所有可能的組合分析形式：使用微生物組和宏基因組數(shù)據(jù)作為輸入，通過MIMOSA2可選構(gòu)建兩個模型，分別為：AGORA和EMBL；本次我為大家演示擴增子數(shù)據(jù)和代謝組數(shù)據(jù)如何構(gòu)建模型并進行后續(xù)分析。從圖中我們可以看到基于Greengenes數(shù)據(jù)庫和SILVA數(shù)據(jù)庫的物種注釋結(jié)果，還有目前最為流行的非聚類方式得到的ASV(amplicon sequence variant)都可以用于分析。這里我們可以看到需要兩個數(shù)據(jù)庫：AGORAh和EMBL數(shù)據(jù)庫。所以首先我們來下載數(shù)據(jù)庫。

下面我們下載數(shù)據(jù)庫，基于非聚類方式ASV的兩種模型數(shù)據(jù)庫和基于傳統(tǒng)聚類(OTU)的兩種模型數(shù)據(jù)庫，注意這里基于OTU的兩個模型適用于greengene和Silva數(shù)據(jù)庫注釋結(jié)果。download_reference_data為數(shù)據(jù)庫下載命令，網(wǎng)站 https://borenstein-lab.github.io/MIMOSA2shiny/downloads.html# 可以手動下載數(shù)據(jù)庫，但是好像我不能下載，顯示拒絕訪問。大家還是用下面命令好些。數(shù)據(jù)庫默認下載到當前文件夾中，所以如果后使用多的話指定一個固定文件夾。四個數(shù)據(jù)庫下載下來大小一共為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”)

載入依賴關(guān)系

```{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)))

下載參考數(shù)據(jù)庫

```{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”)

#基于參考數(shù)據(jù)庫方法

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為模擬的物質(zhì)，M濃度的變化和物種1-3之間的關(guān)系如下圖所示。

B. A-F為不同樣本，物種1和2對于M物質(zhì)具有相同的變化，這表明了這兩種物質(zhì)促進M物質(zhì)形成，物種3消耗M物質(zhì)。比較群落水平的代謝潛力(Compare total community-level metabolic potential, CMP)是最重要的指標，代表每個物種對該物質(zhì)的影響得分。

C. 表示CMP進行回歸并評估對給物質(zhì)影響的顯著性。這里給出R為評估的全部物種對該物質(zhì)的解釋度。

D. 給出了三個物種影響程度，這里最大的`Taxon 3`也就是所謂的支配物種，對M物質(zhì)具有支配作用。

### 運行計算

這里一共有四種模式的運算：

#### 模式1：基于Greengenes OTU的結(jié)果和EMBL模型

第一種模式：使用基于聚類的結(jié)果，也就是OTU進行分析，我們在這里選擇RefSeq/EMBL_GEMs genomes and models模型，其對應(yīng)的數(shù)據(jù)庫再之前下載得到為OTU_EMBL，這里我們指定數(shù)據(jù)庫的相對路徑，并指定到其中的data中。

參數(shù)文件：config_example1.txt，方便理解查查看，這也是我們唯一的輸入，因為全部的參數(shù)都在這里指定好了。這里值得注意的參數(shù)是vsearch_path，我們需要下載vsearch并安裝，這里的安裝無論是再linux還是win下直接賦予可執(zhí)行權(quán)限并放入環(huán)境變量中即可。本例子我運行在win環(huán)境，這里必須寫全稱：vsearch.exe，在linux下只需要寫vsearch就好了。

參數(shù)文件(復制保存為txt即可調(diào)用)：

```{R eval=FALSE, include=FALSE}
file1 test_gg.txt #gg參考數(shù)據(jù)庫聚類結(jié)果
file2 test_mets.txt#代謝物質(zhì)結(jié)果，使用KEGG編號
ref_choices RefSeq/EMBL_GEMs genomes and models#使用模型
file1_type Greengenes 13_5 or 13_8 OTUs#定義輸入微生物數(shù)據(jù)類型
simThreshold 0.99#相似度
data_prefix ./database/OTU_EMBL/data/#定義參考數(shù)據(jù)庫位置
vsearch_path vsearch.exe#定義vsearch位置，此時在win下進行，所以指定為exe
rank_based T
metType KEGG Compound IDs#指定代謝產(chǎn)物ID類型
signifThreshold 1

運行第一種模式

# 下面我們選擇一個例子進行運行，全部輸入文件和參考數(shù)據(jù)庫都備注再這里：config_example1.txt
#---運行-首先測試第一種模式：就gg數(shù)據(jù)庫聚類結(jié)果使用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模型

第二種模式:使用基于聚類的結(jié)果，也就是OTU進行分析，我們在這里選擇AGORA genomes and models模型，其對應(yīng)的數(shù)據(jù)庫再之前下載得到為OTU_AGORA，這里我們指定干數(shù)據(jù)庫的相對路徑，并指定到其中的data中，

參數(shù)文件：./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

```{R}
#---運行=第二個測試基于OTU結(jié)果的RefSeq/EMBL_GEMs genomes and models模型分析
mimosa_results = run_mimosa2("./test4.txt", make_plots = T, save_plots = T)

模式3：基于ASV和EMBL模型

第三種模式: 使用基于非聚類的結(jié)果，也就是ASV進行分析，file1的文件類型需要更改為Sequence variants (ASVs)，我們在這里選擇EMBL模型，其對應(yīng)的數(shù)據(jù)庫再之前下載得到為ASV_EMBL，這里我們指定干數(shù)據(jù)庫的相對路徑，并指定到其中的data中，
參數(shù)文件：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

```{R}
#---運行=第二個測試基于非聚類結(jié)果
mimosa_results = run_mimosa2("./test3.txt", make_plots = T, save_plots = T)

模式4：基于ASV和AGORA模型

第四種模式：使用基于非聚類的結(jié)果，也就是ASV進行分析，file1的文件類型需要更改為Sequence variants (ASVs)，我們在這里選擇AGORA genomes and models模型，其對應(yīng)的數(shù)據(jù)庫再之前下載得到為ASV_AGORA，這里我們指定干數(shù)據(jù)庫的相對路徑，并指定到其中的data中，
參數(shù)文件：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

這里再運行出現(xiàn)錯誤，文件夾作者指定的和數(shù)據(jù)庫中的對應(yīng)不上，我修改了作者的兩個函數(shù)，進行運行：

```{R}
#---運行=第二個測試基于非聚類結(jié)果的AGORA genomes and models模型分析
library(RColorBrewer)#調(diào)色板調(diào)用包
library("cowplot")
library(data.table)
# 主要是build_metabolic_model函數(shù)指定文件夾錯誤，所以我修改了該函數(shù)兩處文件夾的位置
source("./build_metabolic_model-wt.R")
#載入修改后的主函數(shù)run_mimosa2-wt
source("./run_mimosa2-wt.R")
#當我調(diào)整之后，進行分析已經(jīng)沒有問題了
mimosa_results = run_mimosa2("./test2.txt", make_plots = T, save_plots = T)

結(jié)果說明

表格1: 微生物對代謝物的貢獻表格

最重要的一列：varshare，代表了該模型匯總為微生物可以解釋代謝物多少差異。僅僅包括貢獻p值小于0.1代謝物。也就是說這些代謝物受到那些微生物影響，或者那些微生物影響了這些物質(zhì)變化。找出對物質(zhì)變化影響比較大的微生物，作為后續(xù)研究方向。

表格2：原始物種表格匹配到數(shù)據(jù)庫后的新的物種表格

我們無論是基于ASV表格還是參考數(shù)據(jù)庫聚類的OTU表，那種模型，都會重新匹配模擬到新的物種，這個文件展示相關(guān)信息(我們可以看到雖然參加分析的OTU表樣品數(shù)量很多，但是僅展示了和代謝組共有的樣本)

表格3：CMP值矩陣

CMP值矩陣，展示了每個樣本中每種物種對每個物質(zhì)的CMP評估值

表格4：模型統(tǒng)計摘要

模型統(tǒng)計總體簡要統(tǒng)計：參與分析的樣本數(shù)量，物種數(shù)量，參與模型的物種數(shù)量，代謝物質(zhì)數(shù)量，參與模型的代謝物質(zhì)數(shù)量，計算得到CMP值的代謝物質(zhì)數(shù)量，匹配模型的代謝物質(zhì)數(shù)量，顯著的代謝物質(zhì)數(shù)量···

表格5：輸入?yún)?shù)文件，這里包含了我么的輸入文件，模型選擇，數(shù)據(jù)庫指定等信息

圖片文件展示

這里兩種物質(zhì)，對應(yīng)兩個結(jié)果。第一列圖片展示的是多個樣本涉及特定代謝物質(zhì)的CMP值加和做的回歸，R值展示為加粗數(shù)字，數(shù)值越大代表微生物群落的預測能力越準確。第二列是對應(yīng)的微生物的貢獻柱狀圖，柱狀圖右邊的是促進該物質(zhì)生成，左邊的是消耗該物質(zhì)。

手動輸出選擇：以上兩種方式我們得到的默認圖表如果是合并圖例，如果想單獨展示，這里結(jié)果中會包含圖片文件我們可以手動輸出，我在這里提供手動輸出的文件：

#手動輸出CMP擬合代謝物的圖表
p = mimosa_results$CMPplots[[1]]
p+ theme_bw()

#手動導出微生物對代謝物的貢獻圖
p = mimosa_results$metContribPlots[[1]]
p+ theme_bw()

輸出網(wǎng)絡(luò)：

這里我們可以看到基于微生物和代謝物之間建立了一個網(wǎng)絡(luò)，這是一個有向網(wǎng)絡(luò)，varshare值在這里類似我們傳統(tǒng)網(wǎng)絡(luò)中的R值。構(gòu)建網(wǎng)絡(luò)的信息位于varshare值矩陣，varshare值作為邊的權(quán)重文件，我們可以添加對微生物群落對代謝物預測能力的R值。

首先我們來構(gòu)造邊和節(jié)點文件，然后輸出，可以選擇的其他軟件進行網(wǎng)絡(luò)繪制：

# 構(gòu)造邊文件
networ = as.data.frame(mimosa_results$varShares)
# 構(gòu)造一個向量用于承接正負關(guān)系
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)
#構(gòu)造節(jié)點文件，代謝物和微生物共同組成，對分泌物的預測預測效果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 中網(wǎng)絡(luò)可視化(可選)

#構(gòu)造網(wǎng)絡(luò)文件
G #轉(zhuǎn)化為鄰接矩陣，attr：填充權(quán)重
occor.r
#構(gòu)造網(wǎng)絡(luò)文件，network包提供
g # (summary(g))
#做沒有權(quán)重的鄰接矩陣
m #設(shè)置可視化layout
plotcord #添加表頭
colnames(plotcord) = c("X1", "X2")
#添加節(jié)點名稱
plotcord$elements #制作邊文件
edglist edglist = as.data.frame(edglist)
# head(edglist)
#添加邊的權(quán)重
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"))

##這里將邊權(quán)重根據(jù)正負分為兩類
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
#合并之前制作的節(jié)點文件
index = merge(plotcord,node,by = "row.names",all = T)
# dim(index)
# head(index)
index$R = as.numeric(as.character(index$R))
#這里我為了突出代謝物，將全部微生物節(jié)點大小定義為代謝物的十分之一
index$R [is.na(as.numeric(as.character(index$R)))] = min(as.numeric(as.character(index$R)),na.rm = TRUE)/10
plotcord = index
#網(wǎng)絡(luò)出圖
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)

撰文：文濤 南京農(nóng)業(yè)大學

責編：劉永鑫 中科院遺傳發(fā)育所

猜你喜歡

10000+：菌群分析?寶寶與貓狗?梅毒狂想曲 提DNA發(fā)Nature?Cell專刊?腸道指揮大腦

系列教程：微生物組入門 Biostar 微生物組 ?宏基因組

專業(yè)技能：學術(shù)圖表?高分文章?生信寶典 不可或缺的人

一文讀懂：宏基因組 寄生蟲益處 進化樹

必備技能：提問 搜索 ?Endnote

文獻閱讀 熱心腸 SemanticScholar Geenmedical

擴增子分析：圖表解讀 分析流程 統(tǒng)計繪圖

16S功能預測 ? PICRUSt ?FAPROTAX ?Bugbase Tax4Fun

在線工具：16S預測培養(yǎng)基 生信繪圖

科研經(jīng)驗：云筆記 ?云協(xié)作 公眾號

編程模板:?Shell ?R Perl

生物科普:??腸道細菌?人體上的生命?生命大躍進 ?細胞暗戰(zhàn) 人體奧秘 ?

寫在后面

為鼓勵讀者交流、快速解決科研困難，我們建立了“宏基因組”專業(yè)討論群，目前己有國內(nèi)外5000+ 一線科研人員加入。參與討論，獲得專業(yè)解答，歡迎分享此文至朋友圈，并掃碼加主編好友帶你入群，務(wù)必備注“姓名-單位-研究方向-職稱/年級”。PI請明示身份，另有海內(nèi)外微生物相關(guān)PI群供大佬合作交流。技術(shù)問題尋求幫助，首先閱讀《如何優(yōu)雅的提問》學習解決問題思路，仍未解決群內(nèi)討論，問題不私聊，幫助同行。

學習16S擴增子、宏基因組科研思路和分析實戰(zhàn)，關(guān)注“宏基因組”

總結(jié)

以上是生活随笔為你收集整理的中文整合包_MIMOSA2： 基于微生物组和代谢组数据的整合分析的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。